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We examine the spatial extrema (local maxima, minima and saddle points) of the covariant scalars 
(density, Hnbble expansion, spatial curvatnre and eigenvalues of the shear and electric Weyl tensors) 
of the quasi-spherical Szekeres dust models. Sufficient conditions are obtained for the existence of 
distributions of multiple extrema in spatial comoving locations that can be prescribed through initial 
conditions. These distributions evolve without shell crossing singularities at least for ever expanding 
models (with or without cosmological constant) in the full evolution range where the models are 
valid. By considering the local maxima and minima of the density, our results allow for setting 
up elaborated networks of “pancake” shaped evolving cold dark matter over-densities and density 
voids whose spatial distribution and amplitudes can be controlled from initial data compatible with 
standard early Universe initial conditions. We believe that these results have an enormous range of 
potential application by providing a fully relativistic non-perturbative coarse grained modelling of 
cosmic structure at all scales. 


I. INTRODUCTION 

The current era of “Precision Cosmology” has pro¬ 
duced a large amount of high quality observational data 
at all astrophysical and cosmic scales whose theoreti¬ 
cal interpretation requires a robust modelling of self- 
gravitating systems. Conventionally, the large cosmic 
scale dynamics of these sources is examined through lin¬ 
ear perturbations on a ACDM background m, while 
Newtonian gravity (perturbative and non-perturbative 
g], as well as numerical simulations gUS]) is often em¬ 
ployed for self-gravitational systems at smaller galactic 
and galactic cluster scales . 

A non-perturbative approach by means of analytic or 
numerical solutions of Einstein’s equations mzHsi is less 
favoured to analyse observations because (unless power¬ 
ful numerical methods are employed) the high non-linear 
complexity of Einstein’s equations renders realistic mod¬ 
els mathematically untractable. As a consequence, ex¬ 
tremely idealised toy models are used in most cosmo¬ 
logical applications based on fully relativistic and non- 
perturbative methods. The most prominent example are 
the spherically symmetric Lemaitre-Tolman (LT) mod¬ 
els [TUHH] ( see extensive reviews in [31IIH1]), which have 
been used to describe fully relativistic non-perturbative 
evolution of cold dark matter (CDM) sources at all scales: 
from local collapsing structures of galactic scale embed¬ 
ded in an FLRW background (spherical collapse and “top 
hat” models HE]) to Cpc sized cosmic voids in the re¬ 
cent effort to explore the possibility of htting large scale 
cosmological observations without the assumption of a 
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dark energy source or a cosmological constant gsHIS]. 

Even if admitting that LT void models have failed to 
fit large scale observations , non-perturbative gen¬ 

eral relativistic models are still needed and can be useful 
to probe structure formation scenarios and to provide 
theoretical support to cosmological observations within 
the current ACDM paradigm. However, less idealised 
models that are not restricted by spherical symmetry are 
required for this purpose, as the CDM structures we ob¬ 
serve in all scales (from galactic surveys gl] [52] ) is far 
from spherically symmetric. In this context, an improve¬ 
ment on the limitations of LT models is furnished by the 
well known Szekeres solutions [55H55] (see their deriva¬ 
tion and classification in cm IS]), which are still re¬ 
strictive but are endowed with more degrees of freedom. 
While an extensive literature already exists on the usage 
of Szekeres models to address various problems in cos¬ 
mic structure modelling and observations fitting [5BM5] . 
their full potential for theoretical and empirical applica¬ 
tions to Cosmology is still open for future development. 

A specihc issue in which Szekeres solutions (specially 
quasi-spherical models of class I cm) are particularly 
helpful is the modelling of cosmic structure. Previous 
work on this issue [la Ezi EiHm smu] (s ee review in 
0) has addressed the study of various aspects of non¬ 
spherical sources, allowing for a coarse grained descrip¬ 
tion of CDM structures currently observed (see specially 
[niiTiim), which typically consist of spatial distribu¬ 
tions of spheroidal under-dense regions (voids) of typical 
30-50 Mpc size surrounded by a web of elongated fila¬ 
mentary over-dense regions [43] . 

Since a distribution of over-densities and voids can be 
identified with a distribution of spatial maxima and min¬ 
ima of the density, we aim in this work (a complementary 
related work is found in |44| ) at extending and improv- 
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ing this literature by undertaking a comprehensive study 
of the spatial extrema (local maxima, minima and sad¬ 
dle points) of all Szekeres covariant scalars (not just the 
density). 

By working in spatial spherical coordinates, we show 
how the conditions for the existence and location of these 
extrema can be split into interdependent angular and ra¬ 
dial conditions, the latter requiring to be handled numer¬ 
ically. However, by considering Szekeres models compat¬ 
ible with a time preserved Periodic Local Homogeneity 
property based on the existence of a sequence of local ho¬ 
mogeneity 2-spheres (analogous to localised FLRW back¬ 
grounds) , we are able to find rigorous sufficient conditions 
for the existence of distributions of arbitrary numbers of 
these spatial extrema. 

We also examine how spatial extrema can arise by 
means of a “simulated shell crossing” mechanism for 
generic models that do not admit this property. We dis¬ 
cuss extensively the classification of the extrema, the con¬ 
ditions for their avoidance of shell crossings and present 
a procedure to specify their spatial location at all times 
from selected initial conditions. 

These results allow us to set up elaborated networks 
of an arbitrary number of evolving pancake-like over¬ 
densities and density voids, in which the spatial density 
maxima or minima can be placed in prescribed comoving 
spatial locations at all times by realistic initial data. We 
believe that such evolving networks can provide a non¬ 
trivial coarse grained description of observed structures 
in a wide range of astrophysical and cosmological scales 
(specially in the supercluster scale [H]), and as such can 
be a very useful tool in current cosmological research. 

The section by section contents of the paper are de¬ 
scribed as follows: we introduce in section El the Szek¬ 
eres models in spatial spherical coordinates and in terms 
of the q-scalars and their exact fluctuations [37]. The 
Szekeres dipole and its angular extrema are discussed in 
section [hi] (particular cases of Szekeres models from spe¬ 
cial orientation of the dipole are displayed in Table 1). 
The notion of the “LT seed model” is introduced in sec¬ 
tion lYl The conditions that define the location of the 
spatial extrema of all Szekeres scalars are listed and dis¬ 
cussed in section |Vl In section IVTl we introduce Periodic 
Local Homogeneity, showing that this property is a suffi¬ 
cient condition for the existence of an arbitrary number 
of spatial extrema of Szekeres scalars. In section [VTT| we 
show how spatial extrema can also arise in generic models 
not compatible with this property by means of a “simu¬ 
lated shell crossings” mechanism. The classification of all 
spatial extrema as maxima, minima and saddle points is 
found by qualitative arguments in section |VHI| The pos¬ 
sibility of avoiding shell crossing singularities in the time 
evolution of these extrema and possible concavity “inver¬ 
sions” of the extrema (a local maximum evolving into a 
local minimum) are discussed in sections and The 
design of Szekeres models that allow for the description of 
multiple evolving cosmic structures (defined by the spa¬ 
tial maxima and minima of the density) is discussed in 


section XI and a numerical example is examined in sec¬ 
tion XH A summary and our conclusions are given in 
section |XHI| The article contains four appendices: Ap¬ 
pendix ^provides the relation between the metric vari¬ 
ables and the spherical coordinates we have used and the 
standard ones. Appendix |B| lists the regularity conditions 
at the origin and to avoid shell crossings, Appendix 0 
provides a rigorous classification of the spatial extrema 
and Appendix [^ the proof of a necessary condition for 
the existence of spatial extrema of the scalars. 


II. SZEKERES MODELS IN SPATIAL 
SPHERICAL COORDINATES. 


The metric of Szekeres models in spherical coordinates 
is given by I23:Q 
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where ICqo = ICqo{r) is defined in (11), the time depen¬ 
dence is contained in the scale factors: 


vn^ f)n 

a = a{t,r), r = T{t,r) = 1 + —, a'= ^, (4) 

a or 

while the Szekeres dipole W and its magnitude W are 
given by 

W = —X sin 0 cos (j) — Y sin 0 sin (j) — Z cos 0, (5) 

277 TT 

= X'^+ Y^ + Z^ = ^ J J W'^ sm9d0d(l), (6) 

0 0 

where X = X{r),Y = Y{r), Z = Z(r) are arbitrary 
functions satisfying the regularity conditions Ar(0) = 
y(0) = Z(0) = 0 and A'(0) = r'(0) = Z'(0) = 0. 

The main covariant scalars associated with Szekeres 
models are: the density p\ the Hubble scalar "H = 0/3 
with 0 = 'S/aU°‘] the spatial curvature scalar K. = 
(1/6) ^^^77. with the Ricci scalar of the hypersurfaces 


^ By “Szekeres models” we will refer henceforth to quasi-spherical 
models of class I with a dust source (see details on the obten- 
tion and classification of these models in [3S]). The spherical 
coordinates that we use are those defined as a “stereographic” 
projection in [8]. The standard diagonal representation of the 
Szekeres metric and the transformations relating it to 
are given in Appendix [X| 
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arbitrary constant t; the eigenvalues E, ^2 of the shear 
and electric Weyl tensors ('1'2 is the nonzero conformal 
invariant of Petrov type D spacetimes). These scalars are 
expressible in the following concise and elegant form: 


A = A„ + T>^^\ A = p,n,IC, 




Vl/2 = 


(7) 

( 8 ) 


where we have introduced the “q-scalars” Aq and their 
exact fluctuations (for A = p,'H,lC,) defined as 
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where T = -^1 — /C^o &n.d is the determinant of the 
spatial part of the metric 0-(§. The q-scalars are re¬ 
lated to proper volume averages of their corresponding 
scalars with weight factor T [36l [37] . At each 2-sphere 
of constant r the are exact fluctuations around this 
average, which defines an FLRW background as r —)• 00 . 
The q-scalars and their fluctuations are covariant objects 
m and reduce in their linear limit to standard variables 
of cosmological perturbations in the isochronous comov¬ 
ing gauge gSl HI]. 


Evaluating the integral in (101 for each scalar A yields 


the following scaling laws and constraints 
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where the subindex o denotes (and will denote hence¬ 
forth) evaluation at an arbitrary fixed t = tg, thus allow¬ 
ing for the study of the dynamics of the models within 
an initial value framework. Every Szekeres model be¬ 
comes fully specified by initial conditions furnished by 
six free parameters: A, two of {pqo, V-qo, fCqo} plus the 
dipole parameters {X, Y, Z}. The models become deter¬ 
mined either by solving the Friedman equation (12) or by 


^ The q—scalars and their fluctuations have been widely used in 
various applications of LT models I46H49I ISTl . Their properties 
are extensively discussed in these references. In references |88II40| 
the density q—scalar pq is denoted by “pAv” ■ 



FIG. 1. Foliation by non—concentric 2—spheres. Ev¬ 
ery time slice of a Szekeres model is smoothly foliated by 
non-concentric spheres marked by arbitrary constant r with 
surface area dvra^r^. The right hand side figure depicts the 
equatorial projection of the level curves of constant proper 
length (along radial rays) of the 2-spheres of constant r in 
the coordinate diagram in the left hand side figure. Plot¬ 
ted in terms of proper length illustrates how the 2-spheres 
are not orthogonal to the radial rays, which gives rise to the 
non-diagonal metric components in 0-©- The difference 
between the panels of this figure is crucial to understand that 
the apparent rotational symmetry (as in the left panel) seen 
in figures [3|[§ and [^ is a coordinate effect without co¬ 
variant meaning: if these figures were plotted with respect to 
proper radial distances their morphology would look like the 
right panel. 


integrating the evolution equations for the Aq, (see 

[371111]). 

It is very important to emphasise that the spher¬ 
ical coordinates we are using lack a covariant mean¬ 
ing. Hence, “angular directions” indicate locations along 
curves with constant {0,(j)) (radial rays), which are not 
spacelike geodesics as in spherically symmetric space- 
times |50| . The origin worldline [j is not a symmetry 
centre. However, as opposed to the standard dipole co¬ 
ordinates (p, q) (see Appendix 0. the angular coordi¬ 
nates (9,4>) are bounded and thus are very useful to 
specify coordinate locations and for an intuitive under¬ 
standing of the properties of the models, for example: 
it is evident from the metric 0-0 that the surfaces 
of constant t and constant r are 2-spheres, as setting 
dt = dr = 0 yields the metric of a 2-sphere with surface 
areas 47ra^r^. These surfaces constitute a smooth folia¬ 
tion of any time slice by non-concentric 2-spheres (see 


^ We assume henceforth that the LT seed model (section IV i ad¬ 
mits a symmetry centre at r = 0, and thus its associated Szekeres 
models admit an origin worldline whose regularity conditions are 
given in Appendix |B 1| The symmetry centre or origin worldline 
is a physically motivated assumption, but it is not absolutely 
necessary, as regular LT and Szekeres and models exist that ad¬ 
mit either two such worldlines or none (for time slices with the 
topology of a 3—sphere or of a “wormhole” [52]). The results we 
obtain in the present article can be easily extended to this type 
of models (see Appendix D of Ezl). 
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figure]^. This follows from the fact that grr depends on 
the angular coordinates, hence the proper radial length 
along radial rays from r = 0 to points in any 2-sphere 
£ = fg v^grrdr smoothly depends on the angular coordi¬ 
nates of the points. 


III. ANGULAR ORIENTATION: THE 
SZEKERES DIPOLE. 


It follows directly from Q, (|^ and (101 that 

the angular coordinate dependence of the metric and of 
all covariant scalars is entirely contained in the Szek- 
eres dipole W [5]. The spherical coordinates allow us 
to obtain a very precise specification of a unique angular 
coordinate location of the dipole at each 2-sphere of con¬ 
stant r from the extrema of each 2-sphere: the “angular 
extrema” that follow from the condition: 


W.e = = 0, (15) 

which yields two antipodal positions in the (r, 0, (j)) coor¬ 
dinates: 

4 ,- = arctan (^I) , = arccos (^=£_^ 

(16) 

(j)+=TT + (p-, 0+=7r — 0_. (17) 

Since r varies smoothly along the 2-spheres, the solutions 
define in the spherical coordinate system the 
following two curves parametrised by r (see figure : 

B±{r) = [r, 6±{r), (t>±{r)\, (18) 

which cross the origin. The magnitude of W at the curves 
B±{t) is 

W±=W(r,0±(r),</.±(r))=±W(r), (19) 

where the subindex ± denotes (and will denote hence¬ 
forth) evaluation along the curves B±{r) and W is the 
dipole magnitude defined in ([^. Therefore, the angular 
extrema at each surface of constant r necessarily lie in 
their intersection with the curves B±{r) (see hgure [^. 
Since W > 0 along B+{r) and W < 0 along B-{r) we 
have at each 2-sphere of constant r 

Angular maxima of W lie along B+{r), (20) 

Angular minima of W lie along B-{r), (21) 

Saddle point of W at r = 0, (22) 

while allowing for the radial dependence of W the full 
spatial 3-dimensional extrema of W are located in the 
curves B± (r) at some r = such that XX' -(- YY' + 
ZZ' = 0 ^ W = 0. 

Restrictions on the orientation of the dipole in the co¬ 
ordinates {r,6,(j)) correspond to particular cases of Szek- 
eres models that follow from restrictions on the parame¬ 
ters A, y, Z. These particular cases are listed in Table 

1 . 



FIG. 2. Curves of angular extrema. The curves B±{r) 
in the coordinate system (r, 6, (p) defined in ( |18[ ) that follow 
from the angular extrema condition (151 are depicted in blue 
{B+{r)) and red {B-{r)) (colours only appear in the online 
version). These curves determine for every 2-sphere marked 
by arbitrary constant r = ro two points with coordinates 
B+{ro) and B-{ro) where the Szekeres dipole has, respec¬ 
tively, an angular maximum and an angular minimum. These 
two points define from (161-( 171 two antipodal angular direc¬ 
tions [6+(ro), </>+(ro)] andr[0_(ro), (?i-(ro)] associated with the 
angular extrema at each r = ro. The full spatial extrema of 
all Szekeres scalars necessarily lie in points along the curves 
B±{r) at values of r that shift in time. Particular cases of 
Szekeres models occur when the curves B± (r) are constrained 
to a radial ray or a plane in (r, 6, (j>). See Table 1. If the dipole 
parameters are defined as piecewise functions (as in (|71[ )), the 
curves B±{r) become piecewise segments (see figuresand. 


IV. THE LT SEED MODEL. 


We shall use the term “LT seed model” to denote the 
unique spherically symmetric LT model that follows by 
setting W = 0 in the metric 0-(i and in the covariant 
scalars 0-(§. Evidently, every Szekeres model can be 
constructed from its LT seed model by specifying the 
dipole parameters A, A, Z that define W. Therefore, 
Szekeres models inherit some of the properties of their 
LT seed model: the scale factors a, T and the q-scalars 


and 


in (111 are common to both. In fact, the volume integral 
(101 evaluated for LT scalars yields the same q-scalars 
Friedman equation ( [TI| ) and ( [I^ [ISl 137105] . Hence, 


the LT scalars satisfy the same relations 0-§g: 


A(it) — Aq+ 


-(It) 


= -D 


{H) 

(It) > 


4- 


(It) 


- 3 ^(.t): 


(23) 


^ Notice that the density LT fluctuation is exactly the term 

‘^PLT ~ PAv'^ that frequently appears in |38lI4Q| . 
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Dipole oriented along a radial ray: one free parameter. 

Parameter restrictions 

Spherical coordinates 

Rectangular coordinates 

X/0, Y = Z ^0 

hf 

II 

o 

[x, 0, 0] 

Y AO, X = Z ^0 

— 2 ’ ~ 2 ’ 2 ’^ 

[0,1/, 0] 

ZAO, A = y = o 

11 

0 

[ 0 , 0 , 2 ] 

X = aof, Y = bof, Z = Co/ 

tan<()± = ±(&o/ao) 

[aoa;, 6oy, co^j 

two of ao,bo, Co nonzero, / = /(r) 

tan6'± = i-v/fflo + ^ 0/^0 

(general ray) 

Dipole contained in a plane crossing the origin: two free parameters. 

Parameter restrictions 

Spherical coordinates 

Rectangular coordinates 

X,Y AO, Z ^0 


[*, y, 0] 

Z,Y AO, X = 0 

4>± = o,-k 

I0,y,z] 

z,x AO, y = 0 


[x, 0, z] 

a()JC-\-b()Y cqZ = 0 

&0 cot 6± = 

aox + boy + cqz = 0 

ao,bo,co nonzero 

—ao cos 4>± — bo sin ())± 

(general plane) 


TABLE I. Particular cases vs dipole orientation. The angular extrema (161-( 17 \ determine the orientation of the Szekeres 
dipole W given by § at every 2-sphere of constant r. The table displays the dipole parameters X, Y, Z (left column) that 
follow from restricting this orientation (i) to lie along a radial ray and (ii) to be contained in a plane crossing the origin. All 
cases in (i) are equivalent and can be related by rotations in the coordinate system (r, <j)) (the same remark applies to the cases 

(ii)). The cases in (i) seem to correspond to axially symmetric Szekeres models (see section 3.3.1 of [^). The dipole orientation 
is given in spherical coordinates (middle column) and in rectangular coordinates (right column) obtained by the standard 
transformations x = r sinO cos(f>, y = r sin Osin <j), z = rcosO. For these particular cases the spatial extrema of Szekeres scalars 
must lie in their corresponding directions. 


(A) 

with the exact LT fluctuations D;, / related to Szekeres 


fluctuations by: 

_ 


(It) 


r A' 

^(it) “ ^(it) — 
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r-w 
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3r 

(A) 

(It) ■ 


(24) 


Since the LT fluctuations only depend on {t,r), it 
is useful to express the Szekeres scalars A = p, H, JC and 
S, 4'2 exclusively in terms of LT objects and W, either 
as: 

(«) 


/I _ /t I 


TD 


V __ (aL xj, - ^ (‘d 

r-w’ ^ 3 r-w’ 

(25) 

or, alternatively, in terms of and 


A — ^(it)+ 


(a) 

r-w ’ 


E = 


TE 


(It) 


r-w’ 


4*2 = ^ 


(It) 


r-w 

(26) 

with the resulting advantage in both decompositions that 
all the angular dependence is contained in the dipole term 
W. However, the decomposition (251 is more useful, and 
is thus preferable, because the q-scalars satisfy simple 
scaling laws like (11) and thus are easier to manipulate 
than the standard LT scalars We assume hence¬ 

forth that regularity conditions to avoid shell crossings 
(B4|-(B6) hold everywhere, hence 


=0 


^(It) - 


sign of = sign of d[;^J 
(27) 

hold everywhere and thus we can consider the zeroes and 
signs of both fluctuations interchangeably. 


V. LOCATION OF THE SPATIAL EXTREMA 
OF SZEKERES SCALARS. 


The conditions for the existence of a spatial 3- 
dimensional local extremum of A = p, "H, /C at any time 
slice are given by 

A' = 0, A,e=0, A,,^ = 0, (28) 


with similar conditions for E and 4'2 (unless stated oth¬ 
erwise all extrema defined by (28) will be local). We 
examine separately below the angular and radial deriva¬ 
tives. 


A. Angular extrema. 


We obtain from (25) the conditions for the angular 
extrema of A = p, "H, /C at any fixed t: 


’0 T-W 

A r-w 


pggrw,, rA;w,, 
(r-w)2 3(r-w)2 

DjgrW,, rA',W,, „ 

(r-w)2 3(r-w)2 

(29) 


with analogous expressions for E and A! 2 - Evidently we 
have = 0 ^ A^s = A^^ = 0, but the con¬ 

verse is false. However, the condition = 0 defines 
for all (0, (j)) a trivial or degenerate angular extremum of 
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^0- As we explain in section [Vl| this condition leads to 
the Comoving Homogeneity Spheres. Therefore, we have 
the following important result: 

The (non-degenerate) angular extrema of all 
Szekeres scalars A = p,'H, JC, as well as E 
and d> 2 , coincide at all t with the angular 
extrema of the Szekeres dipole W. Since 
= 0 is a necessary (not sufficient) 
condition for (281, the spatial extrema of all 


B. 


Szekeres scalars are necessarily located along 
the curves B±{r). 


Radial location of the spatial extrema: extrema 
of the radial profiles. 


The missing condition in (28) for a spatial extremum of 


A = p, %, ICis, the radial equation A' = 0 for A from (25) 
(and similar equations E' = 0 and = 0 from (261). 


However, the spatial extrema are necessarily located in 
the curves B±{r), hence these radial conditions must be 
evaluated for the “radial profiles” of the scalars evaluated 
in these curves (which for fixed t depend only on r). Asa 
consequence, finding the location of the spatial extrema 
reduces to finding the extrema of the “radial profiles” of 
the scalars given explicitly by: 


A+ = A„ 


TD 


(A) 

(it) 


Tt w 


= ^(It) ± 


WD 


A) 

(It) 


Tt W 


E± = - 


TD 


(«) 

_ i}± 

Tt w’ 


^'2± = 


4^ 

3 TtW’ 


(30) 

(31) 


(A) 

where W and are defined in (|6|) and (|24|), and we 


have used (191 and (25)-(26|. The radial comitions be¬ 
come for A+ and A_ 


A+ = 0 ^ 


A'_=0 ^ 



® The condition dAI = 0 may hold in specific values {t,r) for 
all (see section VII. It implies A = Aq = (from ( |10[ l 

and ( |24[ |), and thus A becomes independent of the angular co¬ 
ordinates. The determinants of the Hessian matrix and all its 
minors (see[^ vanish at points that fulfil = 0, and thus 

this condition marks the same degenerate angular extrema of A 
as a constant function for which all derivatives are trivially zero. 
The equality A = Aq = A^j^^ also holds at coordinate values such 
that W(r, 0,0) = 0, but the derivatives of W do not vanish at 
these values. 


and for and ^ 


^2± 


E(t = 0 ^ 


= 0 ^ 


r(rTW') 

Tt w 

r(r T w') 

Tt w 


n'g + rU” = 0, 


p'n + 1-p'' = 0 . 


(34) 

where we used (24) to eliminate and in terms 


of A'^ and A”. 


The solutions of (32)-(^33jor (34) yield 
the extrema of the radial profiles ([SoIhSTi, wTiich to¬ 


gether with the angular extrema satisfying (29) lead to 
the full 3-dimensional spatial extrema complying with 
(28). Therefore, the location of all spatial extrema of 


Szekeres scalars in the spherical coordinates is given by 


[te±, ’’ei, 0±{re±), <(>±(?'e±)], 


(35) 


where d±, (f>± are given by (@ -([T7| and (teiAei) are 
solutions of the radial conditions (32)-(33) (for extrema 
of A) or ( [M| (for extrema of E and ^7)- We can as¬ 
certain the following points before finding (or discussing 
existence) of solutions of the radial conditions: 

• The regular origin is a spatial extremum for 
every scalar. If standard regularity conditions 
hold (see |B 1|, ([2^ and ([^-(|^ as well as ([M]) 


(and thus (|28[)) 


lold at r = 0 for all t. 


• The spatial extrema of all scalars are not co¬ 
moving. The constraints ((32|)-(33) and (34) de¬ 


pend on time, hence their solutions are (in general) 
different for different values of t and define con¬ 
straints of the form re± = re±{t)- However, to 
simplify the notation we will denote these radial 
coordinates simply by re±. 

• Spatial extrema of different scalars. Evi¬ 
dently, the radial conditions (32)-(33) have differ- 


(32) 

A'^ + rA'^ = 0, 

(33) 


ent solutions for different scalars. 

However, the location and classification of these extrema 
cannot be achieved as long as the main technical difficulty 
remains: finding solution of the radial conditions (32)- 


(33) and (34) for r > 0. Evidently, solving these radial 
conditions, for r > 0 and for generic models and with¬ 
out further assumptions, is practically impossible with¬ 
out numerical methods. However, it is possible to obtain 
sufficient conditions for the existence of such solutions in 
specified comoving radial ranges without actually having 
to solve these constraints. 


VI. SUFFICIENT CONDITIONS FOR THE 
EXISTENCE OF SPATIAL EXTREMA. 

Sufficient conditions for the existence of spatial ex¬ 
trema of Szekeres scalars follows from imposing “local ho¬ 
mogeneity” defined (in a covariant manner) by demand¬ 
ing that the shear and Weyl tensors vanish (which implies 
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S = ^'2 = 0) for all t, but only for selected worldlines or 
surfaces (as opposed to global homogeneity leading to the 
FLRW limit if this condition holds everywhere [71IH1IS7]). 
In fact, such local homogeneity holds at the origin world¬ 
line r = 0 if standard regularity conditions are satisfied 
(see Bl). This type of local homogeneity can also be 
imposed on a discrete set of comoving 2-spheres marked 
by fixed radial comoving coordinate values r > 0: 


Periodic Local Homogeneity (PLH). 

Consider a sequence of n arbitrary nonzero 
increasing radial comoving coordinate values 
and n open intervals between them 


Al = 0 < r < rl, .. ,r^ i < r < r” 


(36) 


A Szekeres model is compatible with PLH in 
the spacetime region bounded by the spher¬ 
ical world-tube (t, r",0,(/>) (see figure if 
at each 2-sphere {rl,9,(j)) the scalars A — 
p, H, K- comply at all t with 


preserved for the whole time evolution [^. We only 
need two of the three scalars to comply with (37) 
and (ISSl) because the fluctuations are related by 


these constraints. 

• They become partitioned by n inhomogeneity shell 
regions (see figure]^ bounded by Co¬ 

moving Homogeneity Spheres. Since condition (37) 


that defines these 2-spheres holds in an asymptotic 
FLRW background [4^, each pair of these con¬ 
tiguous spheres plays the role of time preserved lo¬ 
calised FLRW backgrounds that surround these in¬ 
homogeneity shells where the extrema of the scalars 
are located. 

Since the zeroes and signs of the determine the 
concavity pattern (maxima and minima) of the radial 
profiles of all Szekeres scalars (30)-(31), we have the fol¬ 
lowing result: 


Proposition 1. PLH defined by (37) and 


(38) constitutes a sufficient condition for the 


existence of 2n spatial extrema (not located 
in the origin) of the Szekeres scalars A = 
p, "H, K, and of E, d >2 for all t. 


E, = [ 4 - 2 ]* = 0 

^ [d(^)]. = [d|5]* = [a;]. = o, 

A* = Aqjf = 

(37) 


where the subindex * denotes evaluation 
at the n fixed comoving radii (36). We 
will use the term Comoving Homogene¬ 
ity Spheres to denote the non-concentric 
2-spheres {r\,9,(j)) where local homogene¬ 
ity holds. The origin worldline r = 0 can 
be regarded as the Comoving Homogeneity 
Spheres of zero area. 


Szekeres models admitting PLH are characterised by the 
following properties: 


• They can be specified by initial conditions in which 
any two of the tree initial value functions A^q = 
Pqo, "Hgo comply with 


a;o«) = 0 

^ [Dj;^J]o(r:) = 0 D(^)(r:,0,<^)=O, 

^ Ao{rl)Aqirl) = [A(,)]o(r:), 

^ ^o{rl) = [4/2]o(r:) = 0. 

(38) 


Notice that these initial conditions are sufficient to 
define PLH at all t (pending shell crossings), since 
(|^ and the constraints (13)-(14) and (24) that re¬ 
late the E, 'I '2 and the fluctuations and are 


Besides the extremum at r = 0 (see section 0 , models 
compatible with PLH admit 2n extrema marked by coor¬ 


dinate values (35) with radial coordinates .., and 


distributed in n pairs as follows: 

• the radial coordinates of each pair are (for each 
scalar) solutions of the radial conditions (32|-(33) 
or (34) in the intervals 0 < < 

< r", with the r® given by (36). 


• spatial extrema at rl_^_ 


> 0 are located in the curve 
B+{r) and those at rl_ > 0 in the curve B-{r). 


Proof. As we prove in Appendix]^ PLH defined by (37) 
is a sufficient condition for the existence, at all t, of n ex¬ 
trema of each of the radial profiles A±, E±, [4'2]± in the 
radial ranges 0 < < r*, < r". There¬ 

fore, the radial coordinates of the 2n pairs of spatial 
extrema with r > 0 correspond to the radial coordinates 
rl^± of the extrema with r > 0 of A± and E±, [4'2]± 
given in (30)-(31), as these are the radial profiles of the 


scalars A = p, H, K- and E, 4'2 along the curves H±(r), 
and thus satisfy (29). Since they are solutions of the 
radial constraints (32)-(33) or (34) for all t, they also 


If condition dAI = = 0 holds only for one of the scalars 

A = p, "H, /C, it will not hold for the remaining ones (see exam¬ 
ples for generic LT models in ED)- Since S = ^2 = 0 does not 
hold (though one of these scalars may vanish because of 1^), this 
case does not lead to PLH. Also, = 0 for a single scalar 

does not occur in comoving surfaces and (in general) it only holds 
for restricted ranges of the time evolution. We discuss this case 
in section m 






























(a) 


(b) 


t = u 


t = t, 


FIG. 3. Spacetime diagram of a Szekeres model admit¬ 
ting Peridic Local Homogeneity (PLH). A qualitative 
schematic view of the spacetime evolution of four Comoving 
Homogeneity Spheres represented as circles marked by radial 
coordinates r* that follow from assuming local homogeneity 


conditions (371 that define PLH. The extremum at the origin 


(spatial maximum or minimum) is marked by a thick black 
dot, while four pairs of spatial extrema of a Szekeres scalar are 
marked as blue or red thick dots, depending on their location 
along one of the curves of angular extrema B± (r) depicted by 
blue and red curves (colours only appear in the online ver¬ 
sion). As we show in section [VlII| (see also Appendix [C|) , the 
extrema along B- must be saddle points, while those along B+ 
can be spatial maxima, minima or saddles. The origin world¬ 
line is depicted as a solid black line, whereas the dashed curves 
represent the non-comoving wordlines of the four pairs of spa¬ 
tial extrema with spatial coordinates [rl±, d±(rlj-}, (l>±{rl±)] 
with 6±,(l>± given by (161- 0- The shaded area represents 
one of the comoving shell regions [t, Al,0, ij)] between the Co¬ 
moving Homogeneity Spheres where the extrema evolve for 
all t (pending absence of shell crossings, see section IXI. The 
curves B±{r) become piecewise continuous segments (thus 
providing a more precise angular location of the extrema) for 
dipole parameters given as in ( [tTI (see figures and 1^ 


satisfy (281. Since PLH is preserved by the time evo¬ 


lution, the existence and location of the extrema at the 
intervals A* in the curves B±{r) are also preserved in 
time (pending shell crossings, see section IX). 


Proposition 1 can also be proven by qualitative argu¬ 
ments assisted by figure [i] The following inequalities 
that follow directly from (30)“(31| hold at all t for any 


Szekeres model complying with standard regularity con¬ 
ditions (see Appendixes B 1 and B 21 


+ < < Ag and A(i,) < A_ < Ag if < 0, 

+ > > Ag and > A_ > Ag if > 0, 

('• 



(C) 


■V 




\\ 





0 r 

1 r 

(d) 




FIG. 4. Concavity patterns of radial profiles under the 
assumption of PLH. The panels depict four of the possi¬ 
ble patterns of the concavity of the radial profiles (301 at 


an arbitrary t that follow (qualitatively) from the inequali¬ 
ties (39 )-( |40| | applied to the restrictions that arise from PLH. 
Panels (a) and (b) correspond to keeping the same sign 
in the A*, while (c) and (d) describe the patterns when this 
sign alternates. The prohles of Aq, A+, A_ are depicted 

by the solid black, dashed, blue and red curves (colours only 
appear in the online version), within the three radial inter¬ 
vals a);, a*, a* (shaded rectangles) bounded by the comov¬ 
ing coordinates rl = 0,r);,r*,r*. The extrema of Ag (white 
dots) coincide with the rl, while the extrema of A(it) and 
A± (marked by red and blue thick dots) form a pattern of 
maxima, minima or alternating maxima and minima. Simi¬ 
lar patterns occur for the radial profiles and [4'2]± given 
by These profiles are preserved through the full time 

evolution, pending shell crossings or concavity inversions (see 
sections IX and 0- 


E > 0, > 0 if < 0, 

E < 0, E(iq < 0 if > 0, 

^-2 < 0, < 0 ifD('>)<0, 

^-2 > 0, > 0 if > 0. 


(40) 


( 39 ) 


These inequalities, which also hold for Szekeres models 
compatible with PLH, lead to a qualitative but robust 
inference of the time preserved concavity of the radial 
profiles (30)-([M]) showing various possible sequences of 
maxima and minima along the radial intervals A*. This 
is illustrated in figure]^ 
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• Panels (a) and (b): the fluctuations keep the 
same sign in all A^. There are n saddle points of 
Aq at the r®, the maxima and minima of and 
A- lie inside the intervals and at the r®. 

• Panels (c) and (d): the fluctuations alternate 
signs along the sequence of A^. Maxima and min¬ 
ima of Aq occur at the r®, while maxima and min¬ 
ima of A^, A- lie inside the intervals A^. 

As the figures reveal, the order of the maxima and min¬ 
ima depends on the type of extrema at the origin (or 
the sign of in the first interval Aj). Evidently, any 
combination of the four patterns displayed can also be 
defined by PLH initial conditions ( [M| ). 

PLH allows for the existence at all t of an arbitrary 
number of extrema in the specified intervals A® and 
in specified angles [9±,(j)±], which follow from ([T6|- 0 
through a choice of dipole parameters X, Y, Z. In princi¬ 
ple, even models with an infinite number of extrema can 
be considered for an infinite sequence of r®. However, 
simple dipole configurations (as those found in most of 
the literature [2313113111]) can also be constructed by 
PLH initial conditions (38) that admit a single Comov¬ 
ing Homogeneity Sphere at some r = rl, leading to three 
extrema: one at the origin and the other two located in 
coordinates [re±, 0±(re±), (/>±(re±)] with 0 < re± < r}. 


VII. EXTREMA FROM “SIMULATED SHELL 
CROSSINGS” INDUCED BY THE DIPOLE 
PARAMETERS. 

PLH only provides sufficient (not necessary) conditions 
for the existence of spatial extrema of all Szekeres scalars. 
As discussed in Appendix |D 1[ extrema for monotonic ra¬ 
dial profiles may occur in LT and Szekeres models under 
special assumptions on the profiles. However, as shown in 
the numerical examples of mmiis], spatial extrema of 
all Szekeres scalars may occur, even without the special 
conditions discussed in Appendix |D 1[ in models defined 
by generic initial conditions, and thus admitting com¬ 
pletely arbitrary radial profiles of Aq and of the scalars 


(23) of the LT seed model. 


Generating spatial extrema in generic Szekeres models 
is specially straightforward for extrema along the curve 
B+{r), since (see figurej^ the profiles of A+, E+, ['I' 2 ] + 
grow significantly in a given radial range if the term P—W 
(with W given by in the denominator of (30|-(31) 
becomes sufficiently small for specific choices of dipole 
parameters A, Y, Z, which may occur regardless of the 
choice of initial value functions Aqo- Since P — W = 0 
marks a shell crossing singularity at the curve (see 
B2), and thus the radial profiles diverge 

(which implies that the scalars themselves diverge), then 
a choice of dipole parameters in which P — W is positive 
but close to zero produces a maximum or minimum via 
a large growth (or decay) of these profiles (shown in fig¬ 
ure by creating (in specific radial ranges) conditions 


close to a shell crossing (i.e. regular conditions that sim¬ 
ulate or approximate the behaviour near a shell crossing 
singularity). 

In order to illustrate that a maximum of A+, E+, d(' 2 + 
(and thus a spatial maximum along B+{r)) is easy to 
achieve from a simulated shell crossing, we consider a 
choice of dipole parameters such that P — >V has a 
minimum at some r = rtv > 0 for a fixed t. Let 
r{t, rtv) — VV(Uv) = e > 0 (for sufficiently small e) be the 
minimal value of P — W, hence P'(t,rtv) — W'(r’tv) = 0 
necessarily holds. Under these simple generic assump¬ 
tions the radial conditions (32) and (34) for an extremum 


of either one of the radial profiles A_|_ , S+, ^ 2 + become 
for r « rtv 

A'_^_ «(!-!- 3e)A^(t,rtv) + ftvAq{t, rtv) = 0, 

E+ « -'H'q{t,rtv) - rtv'H"(t,rtv) = 0, 

Att 


%+ ~ Y [p'qit, Uv) + rtvPqit, rtv)] = 0, 


(41) 


and are easily satisfied for a very wide range of generic ini¬ 
tial conditions and radial profiles of Aq in which Aq^O 
holds but A'q and A” have opposite signs, as this com¬ 
bination of signs is very common: Aq grows {A'q > 0) 
with the concavity of a maximum (A" < 0) or Aq decays 
{A'q < 0) with the concavity of a minimum (A" > 0). The 
type of extremum obtained by a minimum of P — W > 0 
at r = rtv depends on the concavity of the profile of Aq, 
which is determined by the sign of and follows from 
the type of extremum at the origin r = 0 that is common 
to A+, E+, 4('2 -i-, but does not depend on the choice of 
dipole parameters because >V(0) = 0. Hence, as illus¬ 
trated in figure if there is a central minimum of A+ 
and ( jdlj ) hold, the extremum at r = rtv is necessarily a 
maximum (and vice versa). 

Under the conditions we have assumed, the radial pro¬ 
files of A_|_, E_|_, 4'2+ in (30)-(31| in a radial range r « rtv 
take the form 




3e 


^tv 'Bq {t, rtv) 

3e ^ 


4' 


2-1- 


47rrtv Pq{t,rtv) 


(42) 


which for e <?; 1 clearly conveys a large growth or decay 
of these profiles depending on the sign of A'q{t,rtv) (or, 
equivalently, the sign of 

An extremum of A_, E_, 4'2_ (and thus a spatial ex¬ 
tremum along B-{r)) can also occur without assuming 
PLH, but the conditions for this are far more stringent 
even if we attempt a simulated shell crossing, since follow¬ 
ing (39)-(|40| (which are valid even without PLH) these 
profiles cannot grow (or decay) more than the LT profiles 
(see panel (b) of figure [^. All this can also be appreci- 
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(a) 



(b) 



(411 change with t, and thus may not hold for the whole 
time evolution (see section 0; (ii) conditions ( [dlj ) must 
be evaluated at all t to verify avoidance of real (not sim¬ 
ulated) shell crossings. All this depends crucially on a 
careful selection of the free parameters and involves de¬ 
tailed numerical work. As a contrast, compatibility with 
PLH effectively guarantees the existence of spatial ex¬ 
trema for all t, though verifying avoidance of shell cross¬ 
ings still remains a technically challenging problem (see 
section IX). 


FIG. 5. Radial extrema from simulated shell cross¬ 
ings. The panels depict the radial profiles A+ (panel (a)) 
and A- (panel (b)) assuming a choice of dipole parameters 
that force a growth of A+ for arbitrary (even monotonic) pro¬ 
files of Ag, A(it) if the term F — W takes small values close 
to zero (a zero marks a shell crossing singularity). We con¬ 
sider the time slice t — to (hence Fq = 1) and use the pa¬ 
rameter choice Z = Y = Q (first entry of Table 1) so that 
yV — X = fcsin^[(5r/18 — l)7r] vanishes at r = 0 and r = 3.6 
with k — 0.5,0.6, 0.75, 0.85 controlling the deviation of 1 — A 
from zero (the vertical dotted line marks r = 1.8 where X is 
maximal). Notice how A+ (panel (a)) exhibits maxima with 
a large growth even if the radial profiles of A, and A(it) are 
both slowly increasing monotonous void radial profiles (we 
used Aq = 0.5 — 0.3/(l -f r^)). 


ated from the expressions equivalent to (41) and (42) 


A' 


1 -h 3e -h 6W(rtv) - 


-I- rtvA"(t,rtv) = 0, 


2rtvW'(rtv) 
2W(rtv) + e 


t’tv) 


A- « Ag(rtv) -h 


^tV Ag (t, Ttv) 

2W(rtv) + e ■ 


(43) 


with analogous expressions for E_, 'I' 2 -. Evidently, as 
shown in panel (b) of figure ( [43| it is much harder to 
satisfy than (41) and the growth of decay of A+ is far 


more restricted than that of A+ in (42). 


Since the extrema that purely follow from simulated 
shell crossings are obtained by (and require) a careful 
choice of dipole parameters, it is very difficult (likely im¬ 
possible) to provide generic sufficient conditions for their 
existence (as we can do by assuming PLH). However, 
PLH depends on the initial value functions Ago common 
to the LT seed model {i.e. conditions (|38|), whereas spa¬ 
tial extrema from simulated shell crossings are completely 
independent of the LT initial conditions, hence generat¬ 
ing such extrema is very useful, for example, in setting 
up Szekeres models constructed from specific LT seed 
models that comply with desired properties (as for ex¬ 
ample LT density profiles compatible with observational 
bounds in [211 ESI |39]). Moreover, the practical imple¬ 
mentation of Szekeres models admitting such extrema 
face the following difficulties: (i) the existence conditions 


VIII. CLASSIFICACTION OF THE SPATIAL 
EXTREMA OF SZEKERES SCALARS. 


Whether we assume PLH or not, the type of the ex¬ 
trema of the radial profile of a scalar (i.e. the max¬ 
ima and minima in figures and does not (necessar¬ 
ily) determine the type of the spatial extremum (max¬ 
ima/minima or saddle point) of the same scalar as a full 
3-dimensional object at each fixed t[^. To properly clas¬ 
sify these spatial extrema (see Appendix [C) we need to 
resort to standard criteria for this purpose based on the 
signs of the determinants of the Hessian matrix and its 
minors evaluated at the extrema location given by (35). 


However, we can also obtain this classification from quali¬ 
tative arguments based on the compatibility between an¬ 
gular maxima/minima and the maxima/minima of the 
radial profiles. 


A. Spatial extrema at the origin. 


It is straightforward to show (see ( |47| ) and 0 that, 
as long as regularity conditions hold (see Appendix |B Ij ), 
all second derivatives of A vanish as r —>■ 0, save for the 
second radial derivative, which in the limit r —>■ 0 for all 
t takes the form: 


\Xr., ~ = 


A" 

^(it) 


- r—0 


(44) 


Since scalars have no angular dependence at r = 0, the 
radial extrema are sufficient to determine the type of spa¬ 
tial extrema. As a consequence of (44), whether we as¬ 


sume PLH or not, the extrema of a Szekeres scalar A at 
the origin r = 0 is of the same type as the extrema of the 
scalar of the seed LT model at the symmetry centre 
r = 0: 


A" 


- r=0 


A" 

^(it) 


- r=0 


< 0 local maximum of (45) 

> 0 local minimum of A. (46) 


^ This fact was overlooked in references 0 138114QH42| which mis¬ 
takenly assume that maxima and minima of the radial profiles 
p± are automatically spatial maxima and minima of the density. 
This assumption is only correct for the extremum at the origin. 
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which implies that either > A (central maximum) 

or < A (central minimum) hold for all points in 

any small neighbourhood around r — 0. If A!' = 0 at 
r = 0, it may be necessary to expand A at r « 0 to 
higher order. 


B. Spatial extrema at r > 0. 


The classification and location of the spatial 3- 
dimensional extrema of A, E, 2 at coordinates (35) with 
Tei > 0 is as follows: 


1. Spatial maxima or minima only occur inside the in¬ 
tervals A* along the curve B+{r), as only in these 
locations the radial and angular maxima and min¬ 
ima coincide. 


2. All extrema along the curve B-(r) are saddle 
points, as we have radial maxima vs angular min¬ 
ima and radial minima vs angular maxima, while 
extrema at r = r® (panels (a) and (b) of figure 
are also saddle points because we have radial max¬ 
ima/minima vs angular degenerate extrema. 

3. Some spatial extrema along the curve B+{r) may 
be saddle points. 

In what follows we prove 1 and 2 above by means of 
qualitative arguments based on the comparison of the 
concavity of the angular and radial extrema. The case 3 
is proven in Appendix [C| 


Classification of the angular extrema. The concav¬ 
ity of the angular extrema of A follows from the signs of 
the angular derivatives along the curves B± (r) 


[AMe± = T 


(he T We) ’ 

(Af+ r,2)[D(^)] 


[A00]e±-T ’ [A«0]e±-o> 


(47) 


where [ ]e± denotes evaluation at and, to simplify the 
notation, Tg and We also denote this evaluation (similar 
expressions arise for the angular derivatives of E and 4'2)- 
The signs of the second angular derivatives above lead to: 


Angular extrema of A along the curve B+ 

minimum if < 0, maximum if > 0, 

degenerate if = 0, (48) 


Classification of the radial extrema. The concavity 
of the radial profiles A± follows from the signs of the 
derivatives A"^. Fortunately, under the assumptions 
of PLH or by generating extrema from simulated shell 
crossings (discussed in previous sections), we can obtain 
the signs of A!^^ directly from the curves displayed in 
figures|^and[^without having to evaluate this derivative. 


Comparison of angular and radial extrema. Com¬ 
bining the information on the signs extrema of A^^± from 
the curves displayed in figures and the information 
from the angular extrema summarised in (48)-(491, we 
obtain for the full 3-dimensional concavity patterns in 
these figures: 


• Assuming PLH: figure 

— Panel (a): In Aj we have > 0, hence 

we have: (i) a radial maximum of A+ (an¬ 
gular maximum) and (ii) a radial maximum 
of A_ (angular minimum). At r = 
we have = 0 and two extrema: ra¬ 

dial minima of A+ and A_ (both degener¬ 
ate angular extrema). The pattern repeats for 
AJ, rj, AJ, rl onwards. 

— Panel (b): We have an analogous pattern of 
extrema as in panel (a), but since < 0 
the types of radial and angular extrema are 
switched: radial and angular minima coincide 
for extrema of A+ inside the intervals A*. For 
all extrema at r = r* we have radial max¬ 
ima and angular degenerate extrema. For the 
extrema of A_ inside the A* we have radial 
minima and angular maxima. 

— Panel (c): All extrema are inside the A* with 
alternating sign. In all odd numbered 
A* radial maxima of A^ coincide with angu¬ 
lar maxima and radial maxima of A_ coincide 
with angular minima. In all even numbered 
A® radial minima of A_|_ coincide with angu¬ 
lar minima and radial minima of A_ coincide 
with angular maxima. 

— Panel (d): As in panel (c), but switching the 
behaviour of odd and even intervals A), (since 
has the opposite alternating pattern of 
signs). 

• Simulated shell crossings without assuming 
PLH: figure 


Angular extrema of A along the curve B- 

maximum if < 0, minimum if > 0, 

degenerate if = 0, (49) 

where the usage of the term “degenerate” for angular 
extrema of Szekeres scalars is explained in section |V A| 


— Since > 0 holds for all the radial range, 
there is a minimum at the origin r = 0 and 
the pattern of spatial extrema is exactly as in 
the first interval A); of panel (c): the radial 
maxima of A^ coincide with an angular max¬ 
ima and radial maxima of A_ (if they occur) 
coincide with an angular minima. 
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— For the case with a maximum at the origin 
(not displayed) we have an analogous situa¬ 
tion as in figure radial and angular minima 
coincide only for A^. 


Comments. The results 1 and 2 emerge readily from all 
the gathered information above. These results also follow 
from the rigorous classification derived in Appendix 
The result 3 follows from the fact (see Appendix!^ that 
the coincidence of radial and angular maxima/minima 
(discussed qualitatively) is a necessary (but not suffi¬ 
cient) condition for a full 3-dimensional spatial max¬ 
ima/minima to occur. However, pending very special 
forms of the free parameters, it is possible to show by 
looking at the behaviour of points in small neighbour¬ 
hoods around the extrema that spatial extrema arise in 
most cases from this necessary condition, with saddle 
points arising when it is not satisfied (along the curve 
B-{r) or at r = r*). 


by 

^(P) ^ Pot)]o ^ ^ [P(S]o ^ rlC^ 

° P 9 O 3pq0 ’ ° Ago 3AgO ’ 

ox(^) _ Q('^)r(p) QC^)rC^) 

ZOq — iigO ^0 ' 

(53) 

(1C’\ 

The constraints on Sq ' depend on the FLRW back¬ 
ground: 



< 1 


if Ago ^ Ao < 0 or ^ fig < 1, 



Open FLRW 

(54) 

j(/C) 

0 ( 1 ) if Ago —> 0 or ^qo —t 0 , 



EdS or ACDM. 

(55) 


Conditions (50)-(52) assure that the initial values of 
Szekeres scalars comply with the expected near FLRW 
conditions at last scattering: po « poj ^0 ~ Hp, Ap « 
Ao and Sq) [^' 2)0 ~ 0. However, notice that (52) does 
not imply W ^ 1, and thus initial Szekeres scalars need 
not be “almost spherically symmetric”. 


IX. SHELL CROSSINGS AND EARLY 
UNIVERSE INITIAL CONDITIONS. 


A. Ever expanding models with A = = 0. 


It is important to verify if the conditions for the ex¬ 
istence of spatial extrema of Szekeres scalars are com¬ 
patible with an acceptable evolution range of the models 
that is free from shell crossing singularities through the 
fulfilment of conditions (B4|-(B6). For general initial 
conditions this is a technically complicated issue that 
must be examined numerically in a case by case basis. 
Hence, in order to obtain useful general conclusions we 
examine below avoidance of shell crossings only under 
the assumption of an asymptotic FLRW background and 
early Universe near homogeneous and near spatially fiat 
initial conditions. 


Assuming as the initial time slice t = to the last scat¬ 
tering surface, the following constraints hold for possible 
FLRW backgrounds: spatially fiat ACDM or Einstein de 
Sitter (EdS), or ever expanding Open FLRW models: 


n 


(m) 

9 O 




1 , n 


(fc) 

qO 


Qfc 
5 2q 


PqO ~ Po, 'HqO « Ho, 


dp) 


< 1 , 




AA«i, 

1-W 


i4 


1 - w 


0 , 


< 1 , 
< 1 , 


qWAIC) 


< 1 , 


(50) 

(51) 


^(fc)|d/C)| 


where -|- — 1 with , H^g ^ defined in 

( |B7| ) and po,iCo,'Ho and are the parameters 

of the FLRW background at the last scattering surface, 
and the dimensionless relative LT fluctuations are given 


Such models are characterised by Ago < 0 ( hype rbolic) 
and Ago = 0 (parabolic). Fulfilment of (|B4[)-( [B^ for all 
t implies the following constraints (which also guarantee 
that p > 0 holds everywhere) that involve only initial 
conditions m 


1 -b slf'> - W > 0 , 1 -b - W > 0 , (56) 

InqoC = - 1 ) + (1 - 4 "^ < 0 , 

(57) 


where we are assuming that 1 — W > 0 holds with W 
defined in ^ and Fq = Hqo{to — Ub) follows from (B 8 ) 
(see analytic form in equation ( 66 ) of [37]). By setting 
= 0 and flg™^ = 1 (which implies Fq = 2/3 from 
(B 8 )) conditions (56) and (57) apply to parabolic models 
(for the latter < 0 must strictly hold for r > 0 ). 


1. Hyperbolic models with a simultaneous Big Bang. 


Since both terms (3/2)Fo — 1 and 1 — Fq in (57) are 
non-negative [IS], setting = 0 (a suppressed decay¬ 
ing mode [46]) implies from (57) the following constraint 
(where we assumed (50)-(52)) 


^(p) _ 

Og — 


4Fn-l 


1 — An 


5 A “gO 


sign 




= —Sign 


On .(/C) 

) ^0 ’ 

(e>) 


(58) 
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which considering that and have opposite signs 
if JCqo < 0, places extra restrictions on the PLH initial 
conditions (38). If the profile of Pqo is like panel (a) of 
figure]^ the profile of ICqo must be also like panel (a), and 
the same goes for profiles like panel (b). If the profile of 
PqO is like panel (c), then the profile of ICqo must be like 
panel (d) and vice versa. 

If besides the general sign conditions on slf \ 
that emerge from ( |58[) we assume early Universe ini¬ 
tial conditions (50 )"(|52p and dipole parameters such that 


1 — W > 0 holds everywhere, we have the following im¬ 
mediate conclusions: 


3. Early time confinement of shell crossings. 

Shell crossings that emerge from the lack of fulfilment 
of (57) {i.e. if > 0) necessarily occur at early times, 
and thus if we assume initial conditions ([50|-(52) these 


singularities can be confined to radiation dominated evo¬ 
lution times much earlier than the last scattering surface, 
where the dust source of a Szekeres model is no longer 
physically acceptable. In order to illustrate this confine¬ 
ment of shell crossings, we assume initial conditions ( |50[ )- 
© and use (0, (0, ([24]) and ( [5^ to examine the 

density in the evolution range t « 


• The first condition in (56) is easily satisfied for any 
FLRW background. 


Fulfilling the second condition in (56) is straight 


forward only for models with an asymptotic Open 
FLRW background (this follows from (51) and (54) 
and (58)). 


The fluctuation Sq ' is necessarily negative (since 
^qo 0 but /Cqo < 0) and is not small (because 
of (55)) at least in the radial asymptotic range. 
Hence, for models with an asymptotic EdS back¬ 
ground the second condition in (56) can be quite 


problematic and its fulfilment requires a careful 
choice of initial parameters. 

The possible shell crossings from lack of fulfilment of the 
second condition in (56) necessarily occur at asymptotic 
late times (see m) , possibly beyond present cosmic time, 
and thus may be irrelevant in usual cosmological applica¬ 
tions (examples of such shell crossings are reported for LT 
models in [Si). The arguments stated above are valid 
whether we assume PLH or not (besides ([50|)- ©)■ 


2. Models with a non-simultaneous Big Bang. 


p,o(l + 4'’^-W) 


P = 


(F - W) 


PgO 


1 + 5 ^^^ - W 


a3/2 (1 + ^(P) _ w) a3/2 - rt'^Jiqo 


PqO 

a3/2 


1 -f - w 


(1 -h - W) a3/2 - 




(60) 


Since the square bracket in the denominator above is 
non-negative (as we assume > 0), a shell crossing 
singularity necessarily occurs as a —>■ 0 at a value 


a = = 


4"’ +1(1 - 

1+4"' - w 


2/3 


> 0 


(61) 


which depends on U\0,4>) but complies with <SC 1 
if conditions (50)-(^ hold. Hence, the shell crossing 
occurs at all spatial coordinates but is entirely confined 
to cosmic times way before the last scattering surface 
where a = Qq = 1. 


The results described above on the first and second 
conditions in (56) remain valid if ^ 0 (irrespective 
of whether we assume PLH or not). Hence, the basic 


difference is the fulfilment of condition (57): 




go^bb 




(p) 




(59) 


which is clearly satisfied by several combinations of radial 
profiles that allow for positive and negative signs of the 
fluctuations (see [5T]). In particular, fulfilment of (59) is 
favoured if < 0 holds, which if we assume PLH 

occurs for profiles of pgo and JCqo given (respectively) by 
the forms of panels (b) and (a) of figure However, 
for profiles with alternating signs of these fluctuations 
displayed in panels (c) and (d) of this figure it may be 
very hard to satisfy (|^. 


4 . Parabolic models. 


We have from (56) and (57) —1 < < 0 

(5(p 1 < 0 for all t [35]. If we assume PLH, then extrema of 
p and T-L (notice that K. — 0) with a complete evolution 
without shell crossings are only possible for radial profiles 
p± and 'H± as in panel (b) of figure For extrema 
obtained from simulated shell crossings on a monotonic 
radial profile, the latter cannot be as in figures and 
but must be clump-like (decreasing with increasing r). 
However, if we assume early Universe initial conditions 
(50)-(5I) at t = to then extrema from void-like profiles 


=> -1 < 


evolving free from shell crossings for all t > to are possible 
because the shell crossings from the violation of (57) can 
be confined to times much before last scattering (section 
IXASj). 
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B. Ever expanding models with A > 0. 


where Hq = 0 and for we have (see equation (74) of 


If A > 0 then ever expanding models follow for initial 
value functions for which the cubic polynomial in the 
denominator of (B7) has no positive roots, which implies 
ICqo < 0 but also allows for some cases with ICqo > 0. 
Since the A term has negligible effect in the early time 
evolution, the behaviour of F — W in the near Big Bang 
range t « is the same as that discussed before for 
ever expanding models with A = 0. However, the A term 
becomes dominant in the asymptotic time range, leading 
from ([B7| to 


r-w« 


3[Kr+ni:'Mr>-iny4<->i 

9 ' 

(62) 

which leads to the following constraints given in terms of 
initial conditions that are implied by (B41-( B 61 


1 + (5^^^ - W > 0, 


1 + - W > 0, 


t'bb < 0 . 


. . 

Since the first and third condition above are identical to 
the first condition in ([5^ and to (57), all the arguments 


concerning these conditions for ever expanding models 
with A = 0 discussed in section IX A| are valid and apply 
directly to ever expanding models with A > 0. However, 
as opposed to the models with A = 0, the possibility 
of late time shell crossings is easily avoided in models 
with A > 0 by selecting early Universe initial conditions, 
since (51)-(52) imply that the second condition in (63) 


can be satisfied for every FLRW background, including 
a ACDM background. As a consequence, if we assume 
initial conditions (50|-(52) with a ACDM background. 


we can avoid shell crossings for all t by restricting initial 
conditions to satisfy < 0 in (59), and if this fails we 


can rest assured (as in the case A = 0) that these singu¬ 
larities are confined to radiation dominated early times 
where the models are no longer valid (section IX A3). 


g'^goCii = ^'^Lqot'bb 




(m) 

90 


[D 


W _ 


^(p) _ y ^ 


3 .(K) 


( 66 ) 

where the analytic form of follows from (57) with Fq 


given by equation (66) of [37] . Evidently, the restrictions 
on the compatibility between PLH initial conditions (38) 


and (64)-(65) are more stringent than in ever expanding 


models. Since initial conditions (38) imply that both 


and have common zeroes with and 5^\ satis¬ 


fying both sign conditions (57) and (66) requires a very 


careful choice of free parameters, more so with the alter¬ 
nating signs of these fluctuations in the radial profiles of 
panels (c) and (d) of figure 


The choice of initial conditions (50|-(52) eliminates 


shell crossings in early times (if = 0 ort(j^ < 0 hold) 
or (if > 0) allow for their confinement well into the ra¬ 
diation dominated times (section IX A 3). However, these 
initial conditions need to be further restricted to satisfy 


the second constraint t' > 0 with given by (|66|). In 


fact, initial conditions (50)-(52) require ' « I also for 


positive spatial curvature, hence the gradient that 


follows from (66) can be very large even if takes its 


linearised form in (59). Therefore, it may be much harder 


to find PLH initial conditions such that shell crossings are 
strictly confined in times that are very close to the col¬ 
lapsing singularity, though collapsing models (or collaps¬ 
ing comoving regions) compatible with PLH may exist 
for which shell crossings are absent in most of the in¬ 
teresting ranges of their evolution. This issue requires 
further study and will be examined in future work. 


C. Collapsing models and regions. 


If A > 0 collapsing models arise for all ca ses with ICqo > 
0 for which the cubic polynomial in (B7) has a positive 


root, if A = 0 then /Cqo > 0 is necessary and sufficient 
for a collapse. Interesting configuration emerge if the 
collapse occurs only in dust layers comprising a region 
around the origin. We examine only the case A = 0, as 
analytic expressions are readily available [371 iS] for the 
conditions to avoid shell crossings: 


I -h Sjf'’ - W > 0, 


Cn > 0 , 


tL < 0 , 


necessary and sufficient, 


(64) 


j(p) _ 


>0, l-t3((5^'’^-h(5^'^^)-W>0, 

only necessary. (65) 


where the second condition in (65) follows from demand¬ 
ing that P — W > 0 holds at the “maximal expansion” 


X. CONCAVITY INVERSIONS OF SZEKERES 
SCALARS. 


We have shown that spatial extrema of all Szkeres 
scalars with specific concavity (maxima, minima or sad¬ 
dle points) arise, either through PLH initial conditions 
(38), or by selecting dipole parameters to induce simu¬ 


lated shell crossings. It is important to verify if the initial 
concavity of these extrema is preserved, or may change, 
during the time evolution (we remark that such concav¬ 
ity inversions occur for the central extremum in generic 
LT models EH). 

We illustrate this concavity change taking as example 
the density fluctuation in the case A = 0 (the discussion 
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is analogous for other scalars and for the case A > 0) 


d(p) 


PqO 

1 

+ - r 




r-w 
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(Fp) 3A/C)\ 

F 2 

PqO 


2^0 ) 

Fq 3 

a3 


r- 

w 


■ 


PqO 


3 (5, 


(p) 


-¥o’"^)^AF-Po)+S', 


(p) 

0 


r-w 


(67) 


where F and Fq are defined i n (|B7[ )-(B8) and we used 
(57). PLH initial conditions (|38[) imply = 

6^ = 0 and =0 for all t at the values r = r\. 

However, the sign of in the intervals between 
the r\ may change with the time evolution. Assum¬ 
ing that regularity conditions (B4)-(B6) hold (absence 
of shell crossings) and a non-negative PqQ, the conditions 
for this follow from the constraint = 0 with 

given by (1^. We have the following cases: 


• If = 0 there are no concavity changes of p in 
the intervals A®. This follows from the fact that 
F/Fq — 2/3 has a definite sign for all t: it is pos¬ 
itive/negative for hyperbolic/elliptic models [46] . 
Hence, the sign of in each A® is preserved by 
the time evolution, as it only depends on the sign 
of ^g'’^ — (3/2)Sg^\ which is preserved for all t (see 
details in [tSjl. 


• If dg^^ and <5g^^ are both nonzero and ^ 0, 
= 0 occurs if 


-Hq 

F.qf) 


(F-Fg) 



( 68 ) 


whose solutions (for 5g^^ ^ §<^ 0 ^^) non- 

comoving values rtv = rtv(t) inside the intervals 
A(, that may introduce concavity inversions in fi¬ 
nite time ranges (these are the “profile inversions 
with turning values” in EH). 


• If <5g^^ = 0 but (5g^^ 7 ^ 0 (with ^ 0), then Dg^^ = 
0 but 7 ^ 0 for t 7 ^ fg: 


r»(p) - PgO F-q F - Fq 

2 ° a^HqoT-W’ 


(69) 


hence, for any given sign of (assumed nonzero 
in the intervals A®), the sign of changes in 
these intervals at t = <o (where F = Fg): it is 
the opposite of the sign of dg^^ for t < to (since 
F — Fg < 0) and is the same as that of (5g^^ for 
t > to (since F — Fg > 0). Hence, a concavity 
inversion must happen at t = to for all A®, but 
the concavity at each interval is mantained for all 
t > to (these are the “profile inversions without 
turning values” in ED)- 


These conditions for concavity inversions of p determine 
(in general) if an initial over-density (local maximum of 
p) can evolve into a density void (local minimum of p). 
Since these conditions involve only quantities of the LT 
seed model ^o^\ F.q/'Hqo, F — Fg), they are the 

same (save for the assumption of PLH) as the conditions 
obtained in EH for LT models with A = 0. The condi¬ 
tions for concavity inversions of the extrema need to be 
separately examined for the other scalars H and /C, for 
models that do not assume PLH and for the case A > 0. 


XI. MODELLING MULTIPLE EVOLVING 
COSMIC STRUCTURES. 

So far we have examined sufficient conditions for the 
existence and the classification of spatial extrema of all 
covariant scalars of Szekeres models (whether compatible 
with PLH or not). We now concentrate on the specific 
case of the spatial extrema of the density, though the 
results we present are directly valid and applicable to 
the remaining scalars. The relevance of density extrema 
follows from the identification: 

• the maximal density value in an “over-density” oc¬ 
curs at a spatial maximum of p, 

• the minimal density value in a “density void” oc¬ 
curs at a spatial minimum of p, 

where by “over-density” and “density void” we mean 
localised regions (“structures”) with significantly higher 
and lower density in comparison with a (suitably defined) 
reference or background density. Evidently, what we call 
“cosmic structure” is roughly a network of such struc¬ 
tures at various scales and levels of coarse graining. In 
particular, we discuss how the spatial location of the den¬ 
sity extrema and their amplitude (relative density growth 
or decay) in these structures can be specified from initial 
conditions. 


A. Location of density spatial extrema when 
assuming PLH. 


1. Radial comoving location. 


PLH constrains the radial coordinates of the spatial 
density extrema in each interval A(. along the curve 
(r) (see figure for the full time evolution (which 
can be free from shell crossings, see section IXI. This 


radial coordinate location is more precise if the values 
that bound the intervals A® are sufficiently close 
(as exemplified in figure]^. While small comoving dis¬ 
tances r® — r®“^ do not (necessarily) imply (at arbitrary 
t) small physical distances in the rest frames (for exam¬ 
ple area distance a (r® — r®“^) or proper radial length), 
a tighter bound on the comoving coordinates of the in¬ 
tervals A® is (without further assumptions) a first step 
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to constrain the radial location of the maximal/minimal 
density in over-densities and voids. 

2. Angular location. 

Since PLH does not restrict the dipole parameters 
X, Y, Z, the latter become free functions to determine 
(through (|I^- 0) the angular location of each spatial 
density extremum at: 0 = cj) = . While 

suitable choices of dipole parameters X, Y, Z can result 
in a curve iB+(r) that covers sufficiently spread angular 
directions in the space (r, 9, tp) (see figures and [§, it 
is more practical to define these parameters in a piece- 
wise manner in order to prescribe a more precise angular 
location at each A*. 

In particular, if the Szekeres dipole W has a single 
free function (see fourth entry in Table 1) then @-0 
lead to curves B±{r) that are radial rays that define two 


constant antipodal angles: (0Oj (po) and (tt — + (po) 

in the coordinate space {r,d,(p). This “axial-like” dipole 
can be obtained if the three dipole parameters are chosen 
as 

X = -cos(pof, Y = -smcpof, Z = -cos6of, 

(70) 

where / = /(r) is an arbitrary non-negative function 
subjected to appropriate boundary conditions, while we 
have introduced the minus sign in order to identify: 

= i0o,(po) and {6_,(p_) = {tt - 9o,tt + cpo), so 
that a choice of {0o,(po) determines the angular location 
of a spatial maxima or minima which are necessarily lo¬ 
cated along S+(r). 

Considering the choice of dipole parameters explained 
above and given n spatial maxima or minima located in 
r = rg_(_ inside radial intervals A*, a convenient form to 
define X, Y, Z is as follows 



r -cos^oi/i, 

Ai 1 

r -sincpoifi, 

Ai 1 

r -cos6)01/1, 

Ai 



1 -COS(p02f2, 

A*, Y=\ 

\ -sin(^02/2, 

Z=\ 

1 -cos 6*02/2, 

A^, 

(71) 

1 

1 (pOn f n •! 

A?, 1 

[ -sintponfn, 

A”, 1 

1 cos 9on /n, 

A?, 



where: 

• the angles Ooi, (poi assign to all r in each radial in¬ 
terval A* (which contain the extrema coordinate 
?’e+) a fixed desired angular direction (see figure]^ 
obtained from 


9+{r)=^o^, (p+(r) = cpoi, r e A(.. (72) 

The curve B+{r) becomes a collection of piecewise 
continuous segments defined in each A® (see figure 
1^. The curve S-(r) is made of such segments lo¬ 
cated in antipodal angular directions (notice that a 
plus sign in 0 simply switches the curve segments 
B+{r) for the curve segments B-[r)). 

• the n non-negative smooth functions fi must sat¬ 
isfy the following boundary conditions to comply 
with smoothness of the scalars (even for piecewise 
continuous curves B± (r)) and to prevent shell cross¬ 
ings: 


/.( 0 ) = /'( 0 ) = /''( 0 )= 0 , 

h^rV) = fM) = = 0 , 

0 < /, < 1 . 

(73) 

Without these conditions the scalars and their first 
derivatives are discontinuous at the r* (as in the 
example in [36)1. 


B. Amplitude. 


While it is possible to control the density growth or 
decay at each over-density or density void in a n interval 
A(. from the choice of initial conditions (38), it is far 


easier to control this growth or decay by the simulated 
shell crossing effect through a suitable choice of dipole 
parameter^ as discussed in section m and illustrated 
by figure Consider a single interval A* in a radial 
profile p+o at the initial time t = to given by (30) and 
following the pattern displayed in panel (a) of figure 
which yields a spatial minimum at r = 0 and a sequence 
of n spatial maxima at each A*. For simplicity we assume 
dipole parameters as in 0 with (poi = 0, 6oi = 7^/2 (first 
entry of Table 1), so that Y = Z = 0 and X = — fi hold 
in A*. A simple convenient form for fi that satisfies the 
boundary conditions (73) is 


fi — hi sin 


{r-rj ^)tt 
ri - rl~'^ , 


W = |A| = /„ (74) 


where 0 < hi < 1, which allows us to use hi to control the 
growth of the density profile p+ in A(. as depicted in panel 
(a) of figurej^ Comparison of this figure with figurej^re- 
veals that a more “narrow” radial range produces a more 
localised simulated shell crossing effect, which allows us 
to control the density growth of decay more effectively. 
Evidently, this process of setting up the amplitude of the 
over-density in figure can be applied for the density ex¬ 
trema in the remaining over-densities at every other A® 
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FIG. 6. Fixing the coordinates of the spatial maxima 
or minima of the density. In Szekeres models compati¬ 
ble with PLH the location of these density extrema in the 
coordinates (r, 6, (f>) can be fixed by: (i) selecting a radial 
interval A* = rl~^ < r < rl, (ii) choosing the dipole parame¬ 
ters so that B+{r) is the thick red line segment in A* (hence 
the angles {Ooi, <j>Qi) specified in (711 are fixed in the interval, 
see further explanation in the text). A network of arbitrary 
number of spatial density maxima or minima (which identify 
over-densities and/or voids) can easily be set up by choosing 
the dipole parameters as in ( |71[ ) so that the curve B+{r) be¬ 
comes a collection of piecewise continuous segments and the 
spatial extrema appear in desired angles in each radial inter¬ 
val A* (see figure Im. This process is also possible for spatial 
maxima and minima of other scalars and for models not com¬ 
patible with PLH involving also a suitable definition of the 
dipole parameters in a piecewise manner. 


(a) 



(b) 



FIG. 7. Amplitude of a density maximum in an over- 
density. The panels depict the radial profiles of p+ and pq for 
the case when PLH is assumed (panel (a)) and when it is not 
assumed (panel (b)), for a maximum located in the interval 
A* = r* < r < r*, with rl = 1.2 and r* = 2.4. We defined the 
dipole parameters as explained in section [XI B| for functions 
fi given by ( 74 1, with k — 0.3 in the intervals A* and A* and 
k = 0.5,0.6,0.75,0.95,0.975 in A* (blue curves). The radial 
profile of Pq in panel (a) complies with the pattern of panel 
(a) of figure]^ and the profile in panel (b) is the same as that 
of figure Notice how a larger value k = 0.975 approaching 
1 in A* yields very large density growth concentrated in this 
radial interval. The dotted vertical line is the maximum of /2 
in A^, which almost coincides with the radial maximum of p+. 
If k = 1.0, this value would mark a shell crossing singularity 
in which diverges. 


XII. NUMERICAL EXAMPLE. 


(and to other Szekeres scalars). In fact, it can be applied 
to every combination of over-densities and voids in the 
A* that follows from the profiles of displayed in figure 
and arbitrary combinations of them. 


Szekeres models contain sufficient degrees of freedom 
to allow for setting up very elaborated and complex net¬ 
works of evolving over-densities and density voids. If we 
assume PLH, the general morphology of these networks 
is also preserved (pending possible shell crossings or con¬ 
cavity inversions discussed in sections IX and 0) by the 
time evolution, as PLH is preserved by the time evolu¬ 
tion (see numerical examples of these evolving structures 

in [13]). 


C. Models not compatible with PLH. 


A. Networks of structures. 


As discussed in section [Vlll spatial extrema of the den¬ 
sity that identify over-densities or voids can be gener¬ 
ated by the simulated shell crossing effect on arbitrary 
radial profiles, leading (for monotonic profiles of the LT 
density) to the two pattern of maxima and minima de¬ 
scribed in section VHI In particular, any regular Szek¬ 
eres model admits dipole parameters defined as in ([tT I, 
hence we also set up networks of over-densities or voids 
in assorted radial and angular locations and can control 
their amplitude by choices of functions fi complying with 
(73). Panel (b) of figure illustrates setting up a den¬ 
sity maximum at a given interval A* when PLH is not 
assumed. 


We provide in this section a simple numerical exam¬ 
ple of a network of structures that is displayed in fig¬ 
ure 1^ Panel (a) of this figure depicts the equatorial 
projection of the level curves of the density contrast 
6 = poIpo — 1 (with po the background density) of an 
initial density configuration at the last scattering surface 
(hence |d| < 10“^). This configuration consists of a large 
spheroidal shaped central void surrounded by five sharply 
defined elongated over-densities obtained from PLH ini¬ 
tial conditions ( [^ with the profile of pqo as in panel (a) 
of figure]^ We selected the dipole parameters as Z = 0 
(0oi = 7r/2 for all i, fifth entry of Table 1), with X, Y 
given as in 0) to determine a unique angular location 
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(a) 


(b) 



FIG. 8. Void with five over—densities. Panel (a) depicts the level curves (in units of 10““^) in the equatorial plane {6 — 'k/2) 
of an initial density contrast 5 = Po/Po ~ 1 at the last scattering surface as a function of y = r/r^ and (j), where po is 
the background density and Ts is a convenient subhorizon scale. We assumed a profile of pqo as in panel (a) of figure for 
dipole parameters chosen as piecewise functions as in (71 1 with 6oi = ■k/2 (so that Z = 0) and all fi given by ( |74| | with 
k = 0.75, 0.8, 0.85,0.9,0.9. The spheroidal density void around the origin is depicted as a blue-cyan area corresponding to 
contrast values —0.3 < 5 < —0.01 (colours only appear in the online version). The void is surrounded by four sharply defined 
over-densities (orange red areas corresponding to 0.01 < S < 1.1), each one located around a density maximum in an interval 
A* in the angular locations (poi = 7r/4, Stt/I, 37r/2, 77r/4, 77r/8. The circles correspond to the Comoving Homogeneity Spheres 
marked by coordinates rl, while the thick black segments are the curves B+{r). The over-densities are filamentary shaped 
in both angular directions (0,(f)), hence their 3-dimensional morphology is pancake like (see figure]^. Panel (b) depicts the 
radial profile of 5 along B+{r) (blue curve) and B-{r) (red curve) together with the profiles of pqo (black dashed curve), the LT 
density contrast (solid black curve) and the zero contrast background value (dotted horizontal line). This type of configurations 
are preserved by the time evolution (see numerical examples in [44]). 


for the spatial density maximum at each over-density in 
the equatorial plane. As shown in panel (b) of figure 
the amplitude of the density contrast for each over¬ 
density is controlled by the radial profile of and by 
the functional form of X and Y (see section VII). 


B. Morphology and evolution of the structures. 

The shape of the central void in panel (a) of figurej^in- 
dicates a spheroidal morphology, which (from qualitative 
arguments) should hold for every structure (over-density 
or void) around a density extremum at the origin. As a 
contrast, the over-densities in panel (a) of figure]^ look 
like elongated thin filaments, a shape that follows from 
the fact that this panel only displays the radial and az¬ 
imuthal {(j) direction) density spread in the equatorial 
projection {6 = 7r/2) of the density contrast, and thus it 
is insufficient to convey the the full 3-dimensional mor¬ 
phology of these over-densities. The missing information 
is provided in figure using as examffie the over-density 
located in the interval A® of figure ^ (the same results 
hold qualitatively for all over-densities not located at 
r = 0). The right hand side diagram of figure shows 
(qualitatively) how the density spreads around the den¬ 
sity maximum along both, the “meridians” (the 9 direc¬ 
tion) and the “parallels” (the (j) direction) of the spherical 


coordinate system. The left hand side plots of this fig¬ 
ure show the angular density profiles (along both angular 
directions) at the fixed radial coordinate of the density 
maximum. Combining the information of these graphs 
and the graph of the radial profile in panel (b) of figure 
we can conclude that the general 3-dimensional mor¬ 
phology of these over-densities is roughly that of “thick 
pancakes”, whose “thickness” depends on the ratio of 
physical distances along the different spatial directions 
in which the density is larger that an given background 
reference value. 


While we obtained the numerical example of figures 
and 1^ by assuming PLH, qualitatively analogous struc¬ 
tures can also be obtained without assuming this prop¬ 
erty, by generating maxima (or minima) of the den¬ 
sity through simulated shell crossings and then selecting 
dipole parameters as explained in section jXlj This proce¬ 
dure to generate a large number of over-densities in the 
equatorial plane was used in [55|. but the dipole param¬ 
eters X and Y (defined as in (711) were not continuous 
in the boundaries separating each piecewise radial range, 
which leads to discontinuous density (a “thin shell” ap¬ 
proximation). As we show in the example depicted in 
figure 1^ this inconvenience can easily be hxed by a suit¬ 
able choice of dipole parameters satisfying appropriate 
boundary conditions (see section]^. 

As a general rule that emerges from figures and 
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FIG. 9. Thick pancake morphology of over—densities. 

The figure illustrates the 3-dimensional “thick pancake” 
shape of the over-densities (with extrema not located in 
r = 0) that arise from Szekeres models. We use as example the 
over-density of figure [^around the density maximum located 
in the interval A® in angular coordinates 605 = 7 r/ 2 , <(>05 = 
Ttt/S. The right hand side diagram shows how this over¬ 
density spreads along both angular directions (j) (red) and 
6 (blue) in the coordinate space {r,6,(j>). The plots on the 
left side show the density contrast as a function of 6 (for 
005 = 77 r/ 8 ) and 0 (for 0 o 5 = 7r/2), with the maximal values 
at the angular coordinates of the density spatial maximum. 
The radial density variation is displayed in panel (b) of figure 
From the angular and radial spread of the density around 
the spatial maximum we can conclude that the over-density 
has roughly a thick paucake morphology. 


the most interesting type of cosmic structure generated 
by Szeleres models (assuming PLH or not) consists of 
a density minimum at r = 0 surrounded by a central 
spheroidal underdense region (density void), which is in 
turn surrounded by an array of curved thick pancake 
shaped over-densities whose density maxima can be lo¬ 
cated in assorted radial and angular directions that can 
be prescribed by initial conditions. Evidently, by using 
the full degrees of parameter freedom described in section 
XIA 2 networks of pancake-like over-densities around 


spheroidal voids can be generated that are more general 
than those of figure with the density maxima in the 
over-densities located all over the 3-dimensional space. 
Also, if we need to construct a mixed network of pancake¬ 
like over-densities and voids, then we can choose for pqo 
and JCgo following different combinations of the radial 
profiles displayed by figure that comply with absence 
of shell crossings at least for the evolution range t > to 
(see section [P^ . 

Figure provides a rough qualitative description of 


the spacetime evolution of initial density configuration 
like the one displayed in figure showing how (under 
the assumption of PLH) each over-density evolves in the 
spacetime shell regions between Comoving Homogeneity 
Spheres. As we showed in section [iXl this evolution can 
be free from shell crossings for a wide variety of early Uni¬ 
verse initial conditions as those used in figure]^ However, 
without assuming PLH it is not evident if a given initial 
concavity (spatial maxima or minima) is preserved for the 
whole evolution, or if concavity inversions may arise (see 
section All this must be examined through a proper 
study of the time evolution of these configurations, which 
requires a separate article (we address these issues in [33] 
and will address them also in future work). 


XIII. SUMMARY AND FINAL DISCUSSION 

We have undertaken a comprehensive unifying study 
of the existence, classification and location of spatial ex¬ 
trema (local maxima, minima and saddle points at hy¬ 
persurfaces of constant t) of the main covariant scalars of 
Szekeres models: A = p,%, K, (density, Hubble expan¬ 
sion and spatial curvature) and the eigenvalues E, '£'2 of 
the shear and electric Weyl tensor. We have expressed 
these scalars in terms of covariant variables |37| (the q- 
scalars and their fluctuations section 0 that reduce in 
their linear limit to standard variables of cosmological 
perturbations [49] . 


A. Radial and angular extrema. 


We have used spatial spherical coordinates that lead 
to a complicated non-diagonal metric, but are very use¬ 
ful for studying the spatial location and the existence 
conditions for spatial extrema of all Szekeres scalars, as 
we can split these conditions (section |v]) into interde¬ 
pendent separate conditions for the angular and radial 
extrema and then look at the full 3-dimensional picture. 
This process lead to the following results: 


• The Szekeres dipole W defines at each 2-sphere of 
constant r a precise angular coordinate location 6 = 
(■'’)) through its angular extrema 

that follow from the condition Wg = = 0. 

For all r these angular extrema define two curves 
(’’) = S±,(j)±] in the coordinate space (r, 0, 0) 

whose functional form at all t is determined by the 
three dipole parameters A, Y, Z (see section HI). 


• The angular extrema of all Szekeres scalars coincide 
with the angular extrema of the Szekeres dipole (see 
. V A and figures and |^. Therefore, finding 


section 


the location of the full 3-dimensional extrema (the 
“spatial” extrema) reduces to solving the radial ex¬ 
trema conditions A(j_ = 0, E(j_ = 0, = 0, where 

7l±, E±, 4>2± are the “radial profiles” obtained by 
evaluating the Szekeres scalars A, E, ^2 along the 
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curves B±{r) (see section VB and figure ffl. In 


general, the radial location of the spatial extrema 
along these curves shifts in time (ie. the worldlines 
of spatial extrema are not comoving). 

• The origin r = 0 (assuming that regularity condi¬ 
tions hold, see Appendix B 1) marks for all times 


the location of a spatial extremum, which (see sec¬ 
tion VIIIAI has the same concavity (local maxi¬ 


mum or minimum) of the extremum of the LT seed 
model at its symmetry centre. 


— a spatial minimum at the origin and an array 
of n spatial maxima. 

— a spatial maximum at the origin and an array 
of n spatial minima.. 

— a spatial maxima or minima at the origin and 
any type of mixed array of n spatial maxima 
and minima. 

• along the curve B- (r) we have necessarily n spatial 
saddle points, 


B. Existence conditions and classification of spatial 
extrema. 


Without solving the radial extrema equations (which 
requires numerical work), we have looked at sufficient 
conditions for the existence of spatial extrema of all 
scalars (section VI and proof in Appendix 0- In par¬ 
ticular, we have proven that arbitrary numbers of spa¬ 
tial extrema for each Szekeres scalar arise by assuming 
compatibility of the models with Periodic Local Ho¬ 
mogeneity (PLH), defined by demanding that the shear 
and Weyl tensors vanish along a sequence of n evolving 
2-spheres (Comoving Homogeneity Spheres, see figure]^. 
PLH is preserved by the time evolution (pending absence 
of shell crossings) and can be specified by suitable (PLH) 
initial conditions that define specific patterns of radial 
profiles (figure]^ which are helpful to understand quali¬ 
tatively the conditions for the existence and classification 
of the spatial extrema. 

We have also examined how spatial extrema of all Szek¬ 
eres scalars can arise without assuming PLH. One pos¬ 
sibility (discussed in section VH and illustrated in figure 

follows by setting up the dipole parameters to induce 
a large growth (or decay) of the scalars in specific radial 
and angular coordinate ranges. We call this procedure 
“generating spatial extrema by simulated shell cross¬ 
ings” , as it involves approximating the behaviour of the 
scalars near a shell crossing singularity, but carefully set¬ 
ting up the free parameters so that the actual singularity 
is avoided. 

A rigorous classification of the spatial extrema 
(whether PLH is assumed or not) is given in Appendix 
[C| However, this classification can also be obtained qual¬ 
itatively (section VHI) by assuming that an extremum is 


a spatial maximum or minimum if its radial and angular 
concavities coincide, while spatial saddles occur whenever 
these concavities do not coincide. For Szekeres models 
compatible with PLH (section [VT| ), the concavity of the 
2n -I- 1 spatial extrema of every Szekeres scalar is given 
as follows: 

• spatial maxima and minima only occur in the origin 
and in each one of the radial intervals A* (between 
Comoving Homogeneity Spheres) along the curve 
H+(r) in the following combinations (see details in 
figure 1^ and in the text of section VI): 


• all extrema located at the Comoving Homogeneity 
Spheres are spatial saddle points. 


For Szekeres models not complying with PLH, whose ex¬ 
trema may follow from simulated shell crossings (see sec¬ 
tion VH I, we have the following combinations that are 
independent of the profiles of the scalars of the LT seed 
model 


• spatial minimum in the origin r = 0 and possible 
arrays of spatial maxima at r > 0 (see figure 

• spatial maximum in the origin r = 0 and possible 
arrays of spatial minima at r > 0. 

These are general conclusions that may admit exceptions, 
since (as we show in Appendix]^ the coincidence of ra¬ 
dial and angular concavities is a necessary but not suf¬ 
ficient condition to determine the concavity of the full 
3-dimensional extrema. 


C. Shell crossings and concavity inversions. 


We examined in section IX the compatibility between 
the conditions to avoid shell crossings (B 2) and the con¬ 
ditions for the existence of spatial extrema. By assuming 
an asymptotic FLRW background and selecting standard 
(near homogeneous near spatially flat) initial conditions 
at the last scattering surface t = to, found the follow¬ 
ing results for models admitting spatial extrema (assum¬ 
ing PLH or not): 


• For ever expanding models (with A = 0 and A > 0) 
an evolution free from shell crossings for all times is 
possible, with some restrictions on the initial condi¬ 
tions for PLH. The exception is hyperbolic models 
with an Einstein de Sitter background in which late 
time shell crossings may occur regardless of the as¬ 
sumption of PLH. 


• Some PLH initial conditions may induce early times 
shell crossings, either with a simultaneous Big Bang 
^bb = 0 or when the condition < 0 is violated. 
These shell crossings can be confined to evolution 
times much before the last scattering, where the 
models are no longer valid (section IX A 3). 
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• Avoidance of shell crossings is much harder to 
achieve for collapsing (elliptic) models or regions 
that admit spatial extrema of Szkeres scalars 
(whether PLH is assumed or not). Shell crossings 
may be eliminated or confined to times before last 
scattering, but it is very difficult at the same time 
to avoid shell crossings in the late collapsing stage 
or near the collapsing singularity. However, rea¬ 
sonable long evolution times without shell crossings 
may be possible. 


Another important technical issue is the possibility of 
concavity “inversions” of the spatial extrema (be. a 
local maximum evolving into a local minimum, section 
H- Some of these inversions occur during a limited time 
range and only for a single scalar (the density but not 
the Hubble expansion or vice versa). Other possible in¬ 
versions occur at a given time slice (see comprehensive 
study for LT models in [ST]). 


D. Modelling cosmic structure. 


We examined the specific case of th e spatial maxima 
and minima of the density (section |XI[ ), which can be as¬ 
sociated with “structures” such as over-densities or den¬ 
sity voids. Whether we assume compatibility with PLH 
or not, Szekeres models allow for the existence of elab¬ 
orated networks of arbitrary numbers of such structures 
around density maxima or minima whose spatial location 
can be specified by initial conditions (section XIA) and 
whose “amplitude” (relative density growth) can be con¬ 
trolled by the choice of dipole parameters following the 
simulated shell crossings procedure (section XIB). 


We have presented a simple numerical graphic exam¬ 
ple (displayed in figures and of a configuration con¬ 
sisting of a central spheroidal void surrounded by over¬ 
densities whose full 3-dimensional morphology is that of 
thick curved pancakes (the combination of central over¬ 
density and pancake voids or any other mixed arrange¬ 
ments are also possible). Since we can prescribe the co¬ 
moving spatial location and amplitudes of the density 
maxima at each over-density in such networks, we do 
obtain a sophisticated toy model that has the potential 
to realise a coarse grained fitting of observed large scale 
cosmic structure imEaiiaiis]. While we have shown 
that such networks are preserved in time under the as¬ 
sumption of PLH, the details of their evolution requires 
separate articles (see [44]). 


It is important to remark that all conclusions we have 
found on the spatial maxima and minima (and saddles) of 
the density apply to the spatial extrema of the remaining 
covariant scalars, which may be useful for a full theoret¬ 
ical understanding of the models and their cosmological 
applications. 


E. Improvements on previous literature. 


There are various elements that intersect the contents 
of our work among the literature on Szekeres models, spe¬ 
cially in references (5] (Sj HSl HZl . However, in sev¬ 

eral key issues our results represent a significant improve¬ 
ment over this existing work. Firstly, all these references 
(including for example [8ll^ IMll?7l[^ 14(114^ 1 have only 
looked at the maxima and minima of the density (over¬ 
densities and density voids), whereas we have considered 
a broad theoretically unifying study of generic spatial ex¬ 
trema (including saddles) of all the covariant scalars of 
the models. We regard this as a specially relevant im¬ 
provement on previous work, as the extrema of the Hub¬ 
ble scalar % or the shear eigenvalue E may be related in a 
non-trivial manner to peculiar velocities and the Hubble 
flow of these structures. Secondly, we have produced a 
clear systematic self-consistent procedure based on rig¬ 
orous sufficient existence conditions and classification of 
an arbitrary number of extrema of all Szekeres scalars, 
as well as on the ability to prescribe their spatial loca¬ 
tion (at all times) of from initial conditions that prevent 
shell crossings at least for the evolution range where the 
models are valid. 


It is important to comment on references [HlinilMlIlO], 
which have also looked at the existence of density ex¬ 
trema and their compatibility with absence of shell cross¬ 
ings. The approach of these references to the density 
extrema is very rudimentary and incomplete: they have 
(specially [10]) only examined the extrema of the radial 
profiles p± (see definition in section V BI and have thus 
assumed (without proof) that these are spatial maxima 
or minima of the full 3-dimensional density (the same 
remarks apply to the intention of identifying structures 
from the “onion model” and the “mock n-body simula¬ 
tion date” of Ensg). However, as we showed in sec¬ 
tion VHI (see rigorous proof in Appendix this is only 
correct for the extrema at the origin, but is completely 
mistaken for the remaining extrema: only the maxima or 
minima of the radial profile (at r > 0 along B+ir)) are 
spatial maxima and minima, while all the extrema of p_ 
(at r > 0 and along B-{r)) are saddles. Here we need to 
point out that the standard coordinates and metric func¬ 
tions used by these references (see Appendix makes 
it much more difficult to obtain a proper classification of 
the spatial extrema of all scalars (including the density). 

Various aspects of our results have been discussed in 
previous articles following a different methodology. Ref¬ 
erences [SHI I3S] also introduced the notion of an LT 
seed model and the dipole parameters are also defined 
as piecewise functions in EZlEni El], but all this has 
been carried on without our systematic and extensive 
approach. Theoretical and observational issues were ex¬ 
amined in Umil EOlEiHSZl on models that describe, 
either an over-density next to a void [261 HZl 133 [39] or 
multiple structures [36| [39] , but the procedure to deter¬ 
mine their spatial location is not systematic and is not 
explained with sufficient methodological rigour. While 
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a nice numerical example of multiple over-densities was 
produced in [36] based on |26l EZ] , the piecewise defined 
dipole parameters were discontinuous, which lead to a 
discontinuous density (as a contrast, we have shown how 
to define these parameters in order to guarantee that all 
scalars are everywhere smooth). The results we have ob¬ 
tained in this article are based on a more solid theoretical 
methodology, and thus greatly enhance and complement 
all this previous literature. 

F. Potential applications and limitations. 

In spite of their limitations (which we discuss further 
ahead), the Szekeres configurations that we have pre¬ 
sented provide a fully relativistic non-perturbative and 
non-linear description of (still) idealised but non-trivial 
evolving CDM structures. The following is a (by no 
means exhaustive) list of potential theoretical and em¬ 
pirical applications for these configurations: 

• Modelling cosmic structure and fitting ob¬ 
servations. The multi-structure configurations 
can be applied to obtain a coarse modelling of ex¬ 
isting large scale structure (from galactic surveys 
EH miiiiiis] or numerical simulations UM), 
thus improving previous work based on spherical 
LT models (whether local structure |9l |55| or Gpc 
sized voids lillSHISl) or on simple Szekeres mod¬ 
els describing an over-density next a density void 
in a crude dipolar array or ‘Swiss cheese” arrays 
[inil271|3H3Sll25|. The networks of over-densities 
and voids that follow from the class I models we 
have examined are also an improvement over the 
plane symmetric “pancake” like structures that 
emerge from Szekeres models of class II |S31[S1] (the 
difference between class I and II models is discussed 
in El)- Our work can also provide theoretical sup¬ 
port to the observational studies of unmiiiisi. 

• Relativistic corrections. By being fully rela¬ 
tivistic and non-perturbative, the multi-structure 
Szekeres configurations are ideally suited to exam¬ 
ine the effects of relativistic and non-linear correc¬ 
tions to Newtonian or perturbative structure mod¬ 
elling or observations fitting nnnsT]. In particular, 
our results can be relevant to complement recent 
work in [411142] . 

• Theoretical issues. By being obtained from an 
exact non-spherical solution of Einstein’s equa¬ 
tions, the multiple structure models are ideal theo¬ 
retical tools to test ideas and concepts that involve 
non-perturbative spacetime inhomogeneity, such as 
cosmic censorship and apparent horizons [55] Unj , 
back-reaction, averaging and the “fitting problem” 
[6nfi64] , as well as the assumption that fulfilling the 
Copernican principle may require large scale homo¬ 
geneity jSS] [50| . The Szekeres multi-structure con¬ 
figurations allow for looking at these issues under 


the assumption of a spacetime inhomogeneity that 
is far less idealised (and more related to the ob¬ 
served structures) than the very constrained spher¬ 
ical inhomogeneity of LT models that is normally 
employed in theoretical studies. 

• The Zeldovich approximation. The pancake 
morphology of over-densities (as those displayed in 
figures 1^ and are reminiscent of structures that 
can be obtained in a Newtonian context through 
the Zeldovich approximation m- In fact, the gen¬ 
eralisation of the latter in General Relativity ISHl- 
|7D| realises the axially symmetric particular case 
of exact Szekeres models of class II (see [SHI [Hij l. 
Since the Zeldovich approximation is an important 
conceptual tool in the study of structure formation 
(even in n-body numerical simulations [HI ED), 
it is important to further explore these theoreti¬ 
cal connections also for the models of class I that 
we have examined. 

• Structure growth and red shift distortion. 

By realising a coarse grained fitting of observed 
large scale structure (galactic clusters, superclus¬ 
ters and void regions) we can achieve a non- 
perturbative and fully relativistic approach to gen¬ 
eralise previous work on these issues [73H75] . 

• Peculiar velocities, Hubble flow and tidal 
forces. Since the locations of the spatial maxima 
and minima of the density are not comoving (see 
figure 3 of [44] and figure [^, we can define and 
study the peculiar velocity field of over-densities 
(along the lines of [53]) with respect to the Hubble 
flow that follows from the level surfaces and spa¬ 
tial extrema of the Hubble scalar H and the shear 
eigenvalue E, while the local field of tidal forces can 
be examined though the level surfaces and extrema 
of the eigenvalue 4 '2 of the electric Weyl tensor. 

• Fitting observations without dark energy. 

Large LT spherical void models used to verify ob¬ 
servations fitting without assuming dark energy or 
cosmological constant have been practically ruled 
out because of the stringent restrictions placed by 
the Kinematic Sunyaev Zeldovich effect [IMS]- 
These restrictions can be re-examined by testing 
this effect on a non-spherical large scale geometry 
based on multiple structures generated by Szekeres 
models (as suggested in [^[5S]). This would im¬ 
prove on previous work using Szkeres models for 
this purpose [28l [29l |32l [39] that were axially sym¬ 
metric or too idealised. 

• Limitations of the models. While Szekeres 
dust configurations are ideal to describe non-trivial 
GDM sources in cosmological scales, they fail to 
incorporate the effects of other forms of cosmic 
matter-energy sources, such as: photon radiation 
(which is not negligible in the early evolution stages 
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near the last scattering time), or internal energy 
in warm dark matter models. Also, Szekeres con¬ 
figurations cannot incorporate the contribution of 
cosmic magnetic fields or gravitational waves, since 
the magnetic Weyl tensor is zero (they are Petrov 
type D solutions). However, all these effects can be 
examined by considering either Szekeres models ad¬ 
mitting sources with non-zero pressure Big, or by 
adding the appropriate perturbative contributions 
associated to these sources to the dust source that 
defines the multiple structures we have constructed 
here. 


from the quadrature in (A3) evaluated at tg (see B8|. 
Then, we re-parametrise the metric functions as: 


a = -, 




zE' 


r=^, W=-, (A4) 


leading to the q-scalars 
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We will undertake the study of these potential applica¬ 
tions and extensions in future publications. 
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Appendix A: Szekeres models in the standard 
representation. 

Szekeres models have been traditionally examined in 
the literature by means of the following metric (see its 
derivation and relation to various parametrisations in 
0 ): 


,2 ,2 ,2 , 2 1 2 ^ 

ds^ = -dt^ + - -^ dz^ + -^(d/ -f dq^), 

1 — K E^ 


E = 


1-b 


p — P 
S 


-b 


q-Q 

s 


(Al) 

(A2) 


where 4) = 4>(t, z), = d^/dz, the functions K, P, Q 

and S > 0 depend only on z and the spatial coordinate 
ranges are z > 0, —oo < p,q < oo. The function 4) is 
determined by the Friedman-like equation 


A2 2M 
4 >^ =-a: -I 

4' 

t tbb = 


Sir 2 


v^d^ 


/o [2M-A:4>-b f^A4-3] 


1/2 ’ 


(A3) 


where M = M{z) and ftb = t^^l,{z) is the Big Bang time 
function. 

In order to relate ( |Al[ ) to the non-diagonal metric (B- 
@ , we choose the z coordinate such that 4’o = 4>(to, 2 :) = 
z holds for an arbitrary t = tg • This choice allows for ex¬ 
pressing tkh as a function of the free parameters M, K, A 


p = P{r) -I- S{r) cot {9/2) cos{(j)), 
q = Q{r) -b S{r) cot (0/2) sin((/), 
z = r, 

(A5) 


which yields the form of W in ([^ with the functions 
X, Y, Z defined as 


l-f. Z.f. (M) 

a re-parametrisation that is justified because the func¬ 
tions S, P, Q only appear in the forms above in the met¬ 
ric and in the covariant scalars 0-0- 

All the results obtained in this article can be fully 
reproduced in the metric representation (011) by apply¬ 
ing the transformations given above. In particular, the 
conditions (38) that define Local Homogeneity Spheres 
(see section VI) consist in demanding that (47r/3)p^g = 
{M/z^y and /C^g = (K/z^)' have common zeroes in a 
sequence of increasing coordinate values z = 0, z)),.., z”. 


Appendix B: Regularity conditions. 
1. Regularity at the origin. 


The following regularity conditions must hold at r = 0 
for all t: 




^(It) = [^\c = 0, 

J C 


D 


M) 


D 


{A) 

(It) 


= 0 , 


PL= [E(tt)], = ['h2L 



= 0 , 


[a]c>0, [a'L = 0, [rL = l, 

Wc = [^1 = = izi = [z'l = 0, 

w. = [W], = 0, 


(Bl) 

(B2) 

(B3) 
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where [ ]c denotes evaluation at r = 0. In particular, we 
have for the dipole parameters X,Y,Z^ 0(r^) for A: > 1 
in any neighborhood r « 0 [5]. These regularity condi¬ 


tions are contained in (37), which justifies describing the 


origin worldline as the Comoving Homogeneity Sphere of 
zero area (see section VI). 


2. Avoidance of shell crossing singularities. 

Since we are only considering models that admit a sin¬ 
gle origin worldline (see Appendix D of m), the regular¬ 
ity conditions to prevent shell crossing singularities are 

[HIST] 


r > 0 

r-w > 0 


(necessary), 

(necessary and sufficient). 


(B4) 


Evidently, condition (B4) restricts the dipole parameters 
X, Y, Z. Since the extrema (angular and radial) of W 
lie in the curves B±{r), then as a consequence of ( |l^ and 
(20)-(21) the following conditions 


r > 0 

r=F >V(r) > 0 


(necessary), 

(necessary and sufficient). 


are equivalent to (B4) but easier to ha ndle because they 
only depend on (t, r). Also, condition (B5| evaluated at 
the initial time slice t = tg (where Tg = I) becomes: 


1 ip yy = I =F Va: 2 -g y2 ^2 > q 


(B6) 


(A3), which takes the following form in terms of q-scalars 

— tbb) = 






Vd dS 




1 1/2 ’ 


(B7) 


where -f - 1, with = 87rp,o/'H^, 


gO 


and = SttA/'H^q, while the Big Bang time 1^,^, is 
given in terms of these initial value functions by 


tbb — to 


F.qO 


1 

F-qO 


Vddd 
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(m) 

go 


-B«m 


B^o’ a3 


-|l/2- 


(B8) 

follows from our choice of 


where Eg = F(l, B^^, B^^^O , 
radial coordinate complying with og = I. 

For the case A = 0 (and some cases with A > 0) the 
quadrature (B7), and thus T and tbb, ca n b e obta ined 
in terms of elementary functions, hence (B4|-(B6) can 


be simple analytic constraints between initial va 


ue func¬ 


tions (see m and in the standard variables in i)- For 


the general case A > 0, (B7) leads to hyper-elliptic inte¬ 
grals and thus (B4)-(B6) needs to be evaluated numer¬ 


ically, with analytic constraints only possible in asymp- 
(B5) totic evolution ranges a « 0 (t « tbO or a I (t tg). 


It is important to remark that (B4) or even its simpli¬ 


fied form (B5) are far more stringent than the condition 
for absence of shell crossings in LT models: F > 0. As 
a consequence, a Szekeres model may exhibit shell cross¬ 
ings even if its LT seed model is free from them. 


which allows us to use absence of shell crossings to con¬ 
strain the dipole parameters. 

Condition (B4) (and its equivalent forms (B5) and 


(B6)) also restrict the free functions of the LT seed model 
because of the presence of T. From its definition in (|^, 
this metric function can be evaluated by implicit radial 
derivative of the quadrature of the Friedman equation 


Appendix C: Rigorous classification of extrema of 
the Scalars. 


At any fixed t a scalar A can be considered a real val¬ 
ued smooth function of {r,6,4>). The classification of its 
extrema as local maxima, minima or saddle points follows 
from the standard “second derivative criterion”: 


Local minimum: det ([HiJ^) > 0, det ([H 2 ]g) > 0, det ([H]g) > 0, 
Local maximum: det ([HiJ^) < 0, det ([H 2 ]g) > 0, det ([HJ^) < 0, 
Saddle point: any other combination of signs. 


(Cl) 

(C2) 

(C3) 


where e denotes evaluation at the extrema: points (re, Oe, 'Pe) that satisfy (28), and H, H 2 , Hi are the Hessian matrix 
and its minors: 


H = 


A^ee 

A,e^ 

A' 


A' 


A^e<f> 

A ,00 A' 

A' 


A 

A" 





, H2 = 

Aae 

A.00 

> 

A ,00 





Hi = [A^ee] ■ 


(C4) 


The second derivatives of A evaluated at an extremum 


located in (28) are the angular derivatives (47) together 


with the mixed radial-angular ones: 

e± 


^TeTWe 


yJAe - A'Fe 
We 


[D(^)]e (X2 + r,2)Z' - i(A2 + i;2)'Ze 


LeTWe 


Wex/AfTi? 
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and the second radial derivative 


tation =pCe is short hand for —C(re+) and C{re-). The 
determinant of the Hessian matrix and its two minors are 


[A" 


\e± 


{-^q)e^e 

TeTWe 


{a± — 


with: a±=- (1 + =p W 


, (C 6 ) 

e 

r T w' 


where /Aq and (to simplify notation) the 

subindex e in the right hand sides denotes evaluation at 
either one of r = re±, while the symbol ± indicates eval¬ 
uation at the curves B±{r): for any quantity C the no- 


det ([Hi] J = T 


(he T We)2 


det([H 2 ]J 

det([H]J 


(X^ + YX ( d [ 5)2 

(TeTWei^ 

(x2 + r/)r2 (d[5)2 

(he T We)2 


(C7) 

> 0, (C 8 ) 

H"±(D[;^))ed/J , 


(C9) 


r. 


with: 'I'p = 


(X'Te - XeY:)^w! + ((X2 + Yl)Z', - i(X2 + Y,^yz,y 


+ Yy)w! 


> 0 , 


(CIO) 


where H" = [A"]^^. Since the determinant of H 2 in (C 8 ) 
is non-negative, we only need to evaluate the signs ot the 
determinants of Hi and H in (C7| and (C9). Hence, the 
second derivative criterion (C1)-(C3) can be stated at 
each extremum B±{re±) as 


local minimum: h±> Q and A± > 0, (Cll) 

local maximum: h± <{) and A± < 0, (C12) 

saddle: any sign combination such that [/i±][A±] < 0, 

(C13) 


with h± and A± given by 


h± = sign [det ([Hi],±)] = T(D[5)e, (CM) 
A± = sign [det ([H],±)] = H" ± (dJ^), vi,,,(C15) 


where we have assumed that Tg > 0 and Tg =F We > 0 
hold from demanding absence of shell crossings. 

Considering the signs of , the nature of the extrema 
above is readily obtained from the following results: (i) 
the combination of signs of the determinants of Hi(H) 
and H(H) given by h± and A± in (C14) and (C15); and 


(ii) the concavity of the radial profiles A± yields the signs 
of the derivative A!' evaluated at the extrema re+, re_ 
(see figure]^. Therefore, the following possibilities read¬ 
ily emerge from (Cll)-(C13) for H depending on the con¬ 


cavity of the radial profiles A± in a given interval A® (see 
figure]^ for reference): 

1. Extremum along H+(r) with radial coordinate r = 

rl+ > 0 : 


(a) If (D|[^^^)e+ < 0, > 0, > 0 then A+ = 

KVI-|(D| 5 )e+ I. We have a local minimum 
if A+ > 0 and a saddle point if A 4 . < 0. 


(b) If (D|[['^j^)e+ > 0, < 0, /i+ <0 then 

A+ = + |(D[;[^^^)e+|. We have a lo¬ 

cal maximum if A+ < 0 and a saddle point if 
A+ > 0. 


2. Extremum along H_(r) with radial coordinate r = 
r®_ > 0 : 


(a) (D|[^j^)e_ < 0, > 0, h_ < 0, hence A_ = 

|A"_| -I- |(D|[['^j^)e_| > 0, therefore only saddle 
points are possible. 

(b) (D|[^j^)e_ > 0, < 0, /i_ > 0, hence A_ = 

— |Wg_| —|(D|[[^^^)e_| < 0 , therefore only saddle 
points are possible. 


These results hold for extrema in any interval A® for all 
combinations of patterns displayed by figure]^ This clas¬ 
sification is easily extended to the extrema of E, '^2 and 
to extrema that arise from simulated shell crossings with¬ 
out assuming PLH (see section VH|. The only extrema 
that are not covered by the criterion (Cl)-(C31 are those 
occurring at the r® where =0 holds (see panels (a) 
and (b) of figure]^. The fact that these extrema are sad¬ 
dles follows from looking at the behaviour of H in a small 
neighbourhood around the points (r®, 0 ±(r®), 0 ±(r®)). 

The rigorous classification of the extrema described 
above coincides with the qualitative one in section |VHI[ 
as expected from the specific sign combinations of h± 
and at B±{r). However, it follows from the form of 
A± in (C15) that even along B-^.{r) same radial and an¬ 
gular concavity and a radial maximum or minimum (sign 
of A"^) are only necessary (but not sufficient) conditions 
for a spatial maximum or minimum of A. 
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Appendix D: Proof of the existence condition of 
extrema of the radial profiles. 

Proposition 1 generalises the following known result for 
LT models (Lemma 3 of [ST] for a single interval Aj = 
0 < r < r^): a zero of at r = r* is sufficient for 
the existence of an extremum of the radial profile of its 
associated covariant LT scalar in Aj. The proof 


follows from the integral form of (24) obtained from the 
integrability of the LT scalar with R = ar: 


D 


(^) 

(It) 


— A(it) Aq — 


4it) dr. 


(Dl) 


Considering the first interval A^ = 0 < r < r* (see figure 
El), we have from ([24|) at any fixed t 


0 = A(it)(t,ri) - Aq(t,rl) 
1 

1 


R^{t,rl) Jo 


^U)JA)R^{t,f)dr, (D2) 


which implies that the integrand must change sign at 
least once in the interval Aj,. Since R > 0 for r > 0, then 
must have (at least) a zero at some value r = rtv 
in this interval. For n intervals A), associated with PLH 
we obtain this result by integrating along each interval 
separately. 

Proposition 1 is proven by generalising Lema 3 of |51j 
to Szekeres scalars and for n intervals A*. For this pur¬ 


pose the LT radial integral (Dl) can be generalised for 
the Szekeres geometry in spherical coordinates as the fol¬ 
lowing 3-dimensional integral 

J j d(j)d9 {Ay^) = J J d(l)d9 j dr [Ay^]' 

(p 6 fp 9 r 

= J J d(l)d9 J dr [A'3;3 + yi(3;3)/]^ 


(p 9 


(D3) 


where the angular variables (0, (f>) take their usual full 
ranges and y satisfies ry'jy = F — W. Considering 
that A'q = (33^'/3^)(A — Aq) follows from ( [l0| ), we have 
A{y ^)' = 33^^3^' Aq-I- A^3^^ = {Aqy^)'. Substitution into 
(|D3|) yields 


d<l>d9 [{A-Aq)y^] =D[5r / I d4>d9 


y^ 


(it )" J J “'1' ““ r _ w 

<p 9 


J J dcl)d9 J dr (A'3^3). 

r 

(D4) 


<t> 9 


where in the second integral we used rA'qjS = (Ap^) — 
A„)r = If we assume PLH and consider the 

9^ (it) - 


first interval AJ, then the left hand side of (D31 vanishes 
identically at r = r*, hence 


0 = J J dct)d9 J dr (A'y^) (D5) 

4 , 9 0 


which implies that A' must have (at least) a zero in this 
interval at some r = rtv that smoothly depends on the 
angular coordinates. In particular, there most be a zero 
of A' at 9±,(j)±. The proof is easily extended to the n 
intervals by integrating along them. 


1. The converse of Proposition 1. 


Notice that the converse of Lema 3 of (51] is false: 
an extremum of A^^) may exist in a given radial range 
r > 0 even without extrema of Aq {i.e. even if we do 
not assume PLH initial conditions (38) so that A^ 7^: 0 
holds everywhere in this range). The proof is as follows. 
We assume that A^ < 0 holds in any given radial range 
r > 0, together with the existence of a zero of at 
some r = rtv > 0, with A'^^^-^{t,r) < 0 for 0 < r < rtv 
and A'^j^^(t,r) > 0 for r > rtv (the proof is analogous 
for the opposite signs). The integral (Dl) for any upper 
integration limit r > rtv can be written as 


D 


{A) 

(it) 





R^df 



A(,j)i?^dr < 0. 


(D6) 

This inequality is not a priori contradictory: it can hold if 
the profile of A(it) sharply decays for r < rtv (for r « rtv) 
and remains nearly constant for r > rtv, hence the first 
integral above can be sufficiently large and negative and 
the second one (which is positive) sufficiently small to 
keep the negative sign of Aq. This result can be easily 
generalised to Szekeres scalars. 
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